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FOREWORD 


This  final  technical  report  entitled,  "A  Deterministic 
Methodology  for  Discriminating  Between  Earthquakes  and  Under- 
ground Nuclear  Explosions,"  is  submitted  by  Systems,  Science 
and  Software  (S 3 ) to  the  Advanced  Research  Projects  Agency 
and  to  the  Air  Force  Office  of  Scientific  Research  (AFOSR) . 
This  report  presents  the  results  of  a research  effort  directed 
toward  the  development  of  an  optimum  multi-discriminant/detec- 
tion  procedure  for  earthquakes  and  explosions.  The  work  was 
performed  under  Contract  Number  F44620-74-C-0063 . Mr. 

William  J.  Best  was  the  AFOSR  technical  contracting  officer. 

Dr.  J.  Theodore  Cherry  was  the  S3  project  manager, 
fc  Thomas  C.  Bache  and  Joseph  F.  Masso  and  Mr.  Bruce 
hason  were  responsible  for  the  development  and  application 
of  the  seismic  ground  motion  prediction  technique.  Dr.  John 
M.  Savino  and  Messrs.  Kenneth  G.  Hamilton  and  David  G. 

Lambert  were  responsible  for  the  analysis  of  the  seismic 
data  against  which  all  theoretical  developments  must 
eventually  be  tested.  Acting  as  consultants  on  the  project 
were  Professors  Charles  B.  Archambeau  of  the  University  of 
Colorado,  David  G.  Harkrider  of  the  California  Institute  of 
Technology  and  Donald  V.  Helmberger  of  the  California 
Institute  of  Technology. 

The  authors  wish  to  extend  their  sincere  appreciation 
to  Ms.  Bernadine  Ludwig  and  Ms.  Darlene  Roddy  for  the  many 
hours  spent  in  the  preparation  of  this  report. 
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I .  INTRODUCTION 


The  fundamental  problem  toward  which  our  research  ef- 
forts were  directed  under  this  project  is  seismic  discrimina- 
tion. Briefly  stated,  the  objective  of  this  program  was  the 
development  of  an  optimum  multi-discriminant/detection 
procedure  for  earthquakes  and  underground  nuclear  explosions. 

Our  approach  to  the  solution  of  the  seismic  discrimina- 
tion problem  was  the  development  of  a deterministic  method- 
ology  for  predicting  teleseismic  ground  motion  from  both 
earthquakes  and  explosions.  This  predictive  capability  pro- 
vides a theoretical  basis  against  which  existing  discriminants 
may  be  tested  and  new  discriminants  designed. 

The  principal  areas  of  research  pursued  during  the 
course  of  this  project  are  the  following: 

1.  Explosion  and  earthquake  source  modeling,  in- 
cluding finite  difference  and  analytic  formula- 
tions. 

2.  Stress  wave  propagation  through  complicated, 
realistic  earth  structures. 

3.  Development  of  state  of  the  art  signal  enhance- 
ment and  identification  techniques. 

4.  Multi-discriminant  design  and  evaluation. 

Each  of  the  sections  of  this  report  is  devoted  to  a 
self-contained  description  of  current  results  relating  to 
one  of  the  four  research  areas  listed  above.  Previous  semi- 
annual technical  reports  written  under  this  contract  (Bache, 
et  al^,  1974,  1975;  Savino,  et  al.,  1975)  have  also  con- 
tained sections  that  are  viewed  as  self-contained  reports 
on  these  and  related  subjects.  We  would  like  to  point  out 
that  this  report  is  not  a summary  of  all  work  carried  out 
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under  this  contract,  but  is  primarily  devoted  to  work  done 
since  the  last  semi-annual  report. 

Much  of  the  work  done  under  this  contract  is  of 
interest  to  the  scientific  community  at  large  and  we  plan 
to  disseminate  the  important  results  by  public*  tion  in  the 
appropriate  scientific  journals.  Two  papers  have  been 
accepted  for  publication  in  the  October  (Bache,  1976)  and 
December  (Bache  and  Harkrider,  1976)  issues  of  the  Bulletin 
of  the  Seismological  Society  of  America.  The  material  of 
Chapters  III,  IV  and  VI  of  this  report  may  be  viewed  as 
drafts  of  papers  to  be  submitted  for  publication. 

The  first  twc  sections  of  the  report  deal  with  the 
three-dimensional  finite  difference  simulation  of  earth- 
quake faulting.  In  Section  II,  a recently  developed  com- 
puter code,  TRES,  is  described.  The  TRES  code  represents 
a finite  difference  simulation  of  stress  wave  propagation 
in  three  space  dimensions.  Special  features  of  this  code 
include  the  capability  to  model  nonhomogeneous  geologic 
environments  and  provisions  for  linking  regions  where 
n^terial  behavior  is  linear  elastic  with  regions  where  the 
behavior  is  nonlinear. 

Section  III  is  devoted  to  a description  of  the  appli- 
cation of  TRES  to  simulate  strike-slip  earthquake  faulting. 

this  three-dimensional  finite  c if  ex  ?nce  code  a numeri- 
cal simulation  of  bilateral  faulting  on  1.0  x 0.6  km  fault 
plane  was  carried  out.  The  basic  mechanism  for  releasing 
the  stored  strain  energy  was  relaxation  of  the  tangential 
stress  at  the  slipping  interface  from  its  welded  value  to 
its  kinetic  friction  value. 

An  important  step  in  our  analysis  of  the  results  of 
the  finite  difference  earthquake  simulation  is  a representa- 
tion of  the  earthquake  generated  elastic  wave  field  in  terms 
of  an  equivalent  elastic  point  source.  This  equivalent 
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elastic  source  representation  is  important  because  it  a] lows 
a linkage  between  the  numerical  codes  and  the  analytical 
wave  propagation  techniques  of  theoretical  seismology. 

Our  primary  interest  in  earthquakes  is  in  their  charac- 
ter as  radiators  of  elastic  waves  to  the  far  field.  An 
advantage  of  the  equivalent  elastic  source  representation 
is  that  it  allows  a convenient  comparison  between  our  finite 
difference  source  and  the  far-field  radiation  characteristics 
of  elastodynamic  earthquake  source  theories.  A number  of 
intriguing  results  are  found  from  a comparison  of  far- 
field  displacement  spectra  from  the  finite  difference  source 
to  those  from  equivalent  earthquakes  modeled  by  the  relaxa- 
tion/volume source  model  of  Archambeau/Minster  (Archambeau , 
1968;  Minster,  1973)  and  simple  dislocation  models  (Haskell, 
1964;  Savage,  1972).  We  find  that  in  almost  all  respects 
(moment,  radiation  pattern,  corner  frequency  behavior)  the 
radiation  from  our  finite  difference  source  more  closely 
resembles  that  from  the  Archambeau/Minster  model.  This  is 
in  spite  of  the  fact  that  in  terms  of  a kinematic  descrip- 
tion of  the  fault  displacements,  the  finite  difference 
earthquake  is  quite  similar  to  a dislocation  model.  In  fact, 
when  we  construct  an  equivalent  dislocation  model  from  the 
outgoing  radiation  field  generated  by  the  finite  difference 
earthquake,  we  find  rather  poor  agreement  between  the 
actual  fault  parameters  (rise  time  and  slip)  and  those  of 
the  equivalent  dislocation  model.  We  attribute  this  dis- 
agreement to  the  fact  that  in  the  finite  difference  calcula- 
tion the  stress  concentrations  that  occur  around  the  fault 
edges  are  relieved  by  plastic  flow.  While  we  do  not  want 
to  place  too  much  significance  on  a single  finite  difference 
calculation,  this  example  does  point  out  a potential  diffi- 
culty with  attaching  physical  significance  to  equivalent 
dislocation  model  parameters. 

In  Section  IV  is  presented  the  results  of  a study  in 
which  we  apply  our  newly  developed  techniques  for  linking 


complex  source  models  (such  as  the  finite  difference  source 
of  Section  III)  with  methods  for  elastic  wave  propagation 
in  realistic  earth  structures.  The  objective  of  the  study 
is  to  compute  synthetic  seismograms  for  a sufficiently  com- 
plex source  model  to  illustrate  the  potential  utility  of 
the  techniques  for  general  seismological  research.  In  doing 
this  we  develop  a scaling  law  relating  teleseisraic  body 
wave  amplitudes  to  the  source  parameters  and  the  elastic 
properties  of  the  material  in  the  fault  region. 

Section  V is  devoted  to  a description  of  the  mode  of 
operation  of  the  Multiple  Arrival  Recognition  System  (MARS) 
computer  code  and  its  application  to  explosion  seismograms . 
The  MARS  code  is  proving  to  be  particularly  effective  for 
detecting,  isolating  and  timing  the  various  seismic  phases 
(P^,  P*,  P^,  , S*,  Sn,  etc.)  that  are  recorded  on  event 

seismograms  in  the  near-regional  and  regional  distance 
ranges.  The  signal  analysis  capabilities  for  narrow-band 
and  polarization  filtering  are  finding  immediate  applica- 
tion to  the  problem  of  depth  determinations  for  shallow 
focus  events.  More  accurate  event  depths  can  be  obtained 
by  incorporating  P and  S wave  arrival  times  in  standard 
hypocentral  location  routines  and  by  the  combined  usage  of 
spectra  of  isolated  body  waves  for  identification  of  spec- 
tral nulls  caused  by  the  interference  of  direct  and  free— 
surface  reflected  phases. 

The  variable  frequency  magnitude  (VFM)  technique, 
embodied  in  the  MARS  code,  has  been  tested  on  both  observed 
and  synthetic  short-period  P wave  seismograms  from  seismic 
events  to  determine  its  effectiveness  as  a discriminant 
between  earthquakes  and  explosions.  In  Section  VI  we  pre- 
sent the  results  of  an  application  of  the  VFM  technique  tc 
a large  population  of  Eurasian  events  recorded  at  LASA  and 
NORSAR  .and  to  North  American  events  recorded  at  the  Yellow- 
knife array  in  Canada.  Obviating  the  requirement  for 
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recording  surface  waves  from  small  explosions  results  in  an 
extended  magnitude  range  over  which  this  discriminant  can 
be  applied.  Based  on  the  event  populations  examined  so  far, 
in  particular  those  occurring  in  Eurasia,  complete  separa- 
tion of  earthquakes  and  explosions  is  achieved  with  the  VFM 
technique  down  to  event  magnitudes,  m^,  in  the  range  4.0  to 
4.5. 

Section  VI  concludes  with  a comparison  of  VFM  esti- 
mates for  theoretical  and  observed  explosion  and  earthquake 
seismograms  for  the  case  where  the  travel  path  effects  are 
comparable.  The  VFM  estimates  for  the  observed  and  synthe- 
tic event  populations  compare  well  in  terms  of  both  absolute 
values  and  population  trends. 


II.  THE  TRES  CODE 


2.1  INTRODUCTION 

The  TRES  code  represents  a finite  difference  Simula— 
tion  of  stress  wave  propagation  in  three  space  dimensions. 

In  the  code  the  conservation  of  linear  momentum,  containing 
the  spatial  derivatives  of  the  individual  stress  components, 
is  written  in  finite  difference  form.  Nonlinear  material 
behavior  is  simulated  by  differencing  the  components  of 
strain  rate  and  incrementing  stress  according  to  an  appro- 
priate constitutive  model.  Explicit  time  stepping  is  used 
to  advance  the  solution  in  time.  The  output  of  the  code 
is  usually  time  histories  of  displacement,  velocity,  stress 
and/or  their  derivatives  at  selected  distances  from  the 
source. 


The  limitations  of  the  code  are  as  follows: 

1.  The  code  is  written  in  three  dimensional 
Cartesian  geometry  and  the  basic  element 
is  a rectangular  parallelepiped. 

2.  All  Lagrangian  rates  of  change  are  replaced 

by  local  rates  of  change.  This  small  deforma- 
tion assumption  linearizes  the  differential 
equation  governing  the  conservation  of 
momentum  and  removes  the  distinction  between 
Eulerian  and  Lagrangian  coordinates. 

3.  Stability  of  the  transient  solution  requires 
that  the  time  step  (At)  be  less  than  the 
time  required  to  propagate  a compressional 
wave  through  an  element  in  the  grid. 

The  code  contains  the  following  special  features : 

1.  Contact  discontinuity  boundary  conditions  are 
used  to  obtain  the  normal  and  tangential 
components  of  stress  over  a specified  interior 
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surface  of  the  computational  grid.  If  the 
tangential  stress  is  then  relaxed  at  a point 
on  the  surface,  a discontinuity  in  tangential 
velocity  results.  This  technique  has  been 
used  to  simulate  the  earthquake  presented  in 
Chapter  III. 

2.  Both  exterior  and  contact  discontinuity  boundary 
conditions  follow  naturally  from  the  interior 
difference  equations  used  for  conservation  of 
momentum.  This  feature  of  the  difference  equa- 
tions also  permits  nonuniform  zoning  to  be  used 
in  the  grid. 

3.  A constitutitve  equation  is  able  to  be  specified 
for  each  element  in  the  grid.  Therefore  a non- 
homogeneous  geologic  environment  is  easily 
modeled.  Furthermore,  the  code  contains  a 
provision  for  linking  regions  where  material 
behavior  is  linear  elastic  with  regions  where 
the  behavior  is  nonlinear.  The  divergence  and 
curl  of  the  displacement  field  may  be  monitored 
over  a spherical  surface  in  the  elastic  region. 
The  equivalent  elastic  source  discussed  in 
Chapter  III  is  obtained  from  these  derivatives 
of  the  displacement  field. 


2. 2 CONSERVATION  OF  LINEAR  MOMENTUM 

If  u,  v,  and  w are  the  x,  y,  and  z components  of 
velocity,  then  the  conservation  of  linear  momentum  may  be 
written: 


du 

dt 
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+ 
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xz 
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(2.2) 


dv 

dt 


iz^  + iiLliZ 

p3x  p9y 


3a 

+ 

p3z 


dw  _ ^axz  ^ayz  ” ax 

dt  p3x  p3y  p 3 z 


(2.3) 


There  will  be  eight  elements  surrounding  an  interior 
nodal  point.  The  centroids  of  these  elements  are  labeled 
A through  H as  shown  in  Figure  2.1.  The  following  incre- 
ments in  x,  y,  and  z give  the  distances  from  the  interior 
nodal  point,  at  which  the  acceleration  is  to  be  computed,  to 
the  corresponding  cenLroid. 


Ax+  A,D,E,H 
Ax”  -*■  B,C,F,G 


Ay+  A,B,E,F 
Ay"  -*■  C,D,G,H 


. + 

Az  -*•  A,  B ,C , D 
Az”  -v  E,F,G,H 


If  Z is  a typical  stress  component  in  the  above  equa- 
tions, then  the  differenced  form  of  the  spatial  derivatives 
used  in  the  code  are  as  follows: 

31  _ ZAB(1  - syHi  - Cz)  vey(i  - iz) 
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Figure  2.1.  The  centroids  of  the  eight  elements  surrounding 
an  interior  nodal  point  are  A,B,C,D,E ,F,G ,H. 
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The  weighting  factors  in  the  above  agnations  are: 


Cx  = 


Ax’ 


Ax+  + Ax 


r,  = 


- +Ay+  rz  = . 

Ay+  + Ay"  Az  + + Az 


Az 


r ' (2.7) 


where  pA  is  the  density  of  the  A element  and 


ZAB  = ZA  " ZB‘ 


(2.8) 


Each  stress  difference  is  weighted  based  on  its  location  rela- 
tlVS  t0  the  interior  nodal  point  at  which  the  acceleration 

rela- 


. - — ,,m*w**  wic  auueier 

IS  to  be  computed.  For  exampie,  the  coordinates  of  I 
trve  to  the  interior  nodal  point  are  (0,  Ay+,  Az+) . “e 


r — ay  , Az  ).  The 
or  - 5 ) (1  - f ) extrapolates  the  £ acceleration 
term  to  the  interior  nodal  point. 

2 • 3 BOUNDARY  CCNDlTTnMc: 

If  a nodal  point  is  located  at  a boundary  then  only 
certarn  elements  surrounding  the  point  will  contain  stress 
It  is  these  real  elements  along  with  the  stress  at  the 
boundary  which  contribute  to  the  acceleration  of  the  nodal 
point.  The  weighting  factors  given  in  Eg.  (2.7)  provide 
* convenient  technioue  for  isolating  the  real  elements. 

For  example,  suppose  that  the  x,  y plane  is  a grid 
boundary  and  the  real  elements  are  in  the  positive  , direc- 
tion, i.e. , A,  B,  c,  D in  Figure  2.1.  If  = 0 in  Egs. 

(2.4)  and  (2.5),  then  only  real  elements  contribute  to  these 
acceleration  terms,  if  E*  is  used  to  denote  the  stress 

the  boundary,  then  the  first  term  in  Eg.  ,2.6,  may  be  written 


£ae(1  ~ 5X)  (l  - ;y)  (EA  - I*)  ,1  - 5X)  (1  . ?y, 


p.Az  + p Az 
A 


PaAz 
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Similar  expressions  involving  I»,  £*,  £*,  pB,  P(;,  Pq  result 
for  the  remaining  terms  in  Eg.  (2.6). 

£ = 1 in  Eqs . (2.4)  and  (2.5),  then  the  real 
elements  are  in  the  negative  z direction,  i.e.,  E,  F,  G, 
and  H in  Figure  2.1.  The  first  term  in  Eq.  (2.6)  is  written 

Eae(1  - 5X>  (1  - 5y>  <*i  - rF)  a - 5X)  (1  - 5y) 

T+ -= ; — = (2.10) 

pAAz  + pEAz  PeAz 

with  similar  modifications  for  the  remaining  terms. 

At  an  exterior  boundary  the  £*  stresses  are  specified 
in  order  to  provide  symmetry  conditions  to  simulate  a free 
surface  or  to  initiate  the  stress  wave.  At  an  interior 
boundary  the  £*  stresses  are  first  obtained  by  assuming 
continuity  of  acceleration  at  a nodal  point  on  the  boundary, 
i.e.  , 


du+ 

du 

dt 

dt  ' 

(2.11a) 

dv+ 

dv 

dt 

dt  ' 

(2.11b) 

dw+ 

dw 

dt 

dt  * 

(2.11c) 

If  the  boundary  is  the  x,y  plane,  then  all  + derivatives 

are  evaluated  with  £Z  = 0 and  all  - derivatives  are  evaluated 
• z 

with  £ =1.  If  we  assume  that 


* 

A 


Ec 


(2.12) 


then  continuity  of  stress  at  the  boundary  along  with  Eqs. 
(2.11)  permit  a unique  solution  for  a*  , a*  , and  a*  at 
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tl-ie  nodal  point.  If  these  stress  components  are  used  to 

move  the  nodal  point,  then  the  acceleration  will  be 

continuous  at  the  boundary.  Relaxation  of  either  c*  or 
^ X z 

°yz  violate  Eqs.  (2.11a)  or  (2.11b)  and  shear  failure 

(slipping)  is  simulated.  Relaxation  of  o*  simulates  ten- 
sion failure. 


2.4  THE  CALCULATION  OF  VELOCITY,  DISPLACEMENT  AND  THE 
TIME  INCREMENT 

the  code,  all  Lagrangian  rates  of  change  are 

replaced  by  local  rates  of  change.  This  assumes  that 
d d d s 

u dx  + v dy  + w dz  is  sma11  coraPared  to  j~.  The  fundamental 

consequence  of  this  assumption  is  that  a given  element  re- 
mains a rectangular  parallelepiped  during  all  states  of 
deformation.  At  a nodal  point,  velocities  are  incremented 
as  follows: 


un+*i  „ un-*l 

vn+%  _ vn-ti 
wn+h  . „n-% 
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3u  A.n 
St 


3v  ..n 
St 


3w  . . n 

^ At 


9 


9 


(2.13a) 

(2.13b) 

(2.13c) 


where  superscripts  indicate  cycle  number, 

Atn  is  the  time  increment  at  cycle  n, 

and  all  acceleration  components  are  evaluated  from  the 
finite  difference  approximation  to  Eqs.  (2.1),  (2.2)  and 
(2.3)  . 

The  corresponding  components  of  displacement,  U,  V, 

W are 


U 


n+1 


Un  + un+!s  itn+!<  , 


(2.14a) 
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vn+1  = vn  + vn+,s  Atn+Js, 

wn+1  = W°  + Wn+Js  Atn+Js  . 

If  ALn  is  the  minimum  dimension  of  all  the  elements  in  the 
grid,  then 

A4-n+35  ALn 

At  = a , (2.15) 

CL 

where 


(2.14b) 

(2.14c) 


« 


I y 


= the  compressional  wave  velocity. 


a is  an  input  constant  s 0.5, 


Then 


At“  - 


At 


n-h 


(2.16) 


2.5  STRAIN  RATE 

After  the  calcula  i in  of  acceleration,  velocity  and 
displacement,  the  code  proceeds  to  calculate  the  strain 
rate  produced  in  the  element.  The  centroids  of  the  eight 
elements  surrounding  an  interior  nodal  point  are  as  shown 
in  Figure  2.1.  The  exterior  points  of  element  A are  labeled 
1 through  8 as  shown  in  Figure  2.2.  The  following  notation 
has  been  adopted. 

u = u - u , x = x -x,  etc.  , 

12  1 2 12  1 2 

, x , y . z 

Ax  = -A^-,  Ay  = Az  = 


25 


The  differenced  form  of  the  strain  rate  components  are  as 
follows : 


L = — 
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(2.17) 


**=  [4xx  ’If1’ 


(2.18) 


* ’ r*  IV, 
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(2.19) 
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(2.20) 


(2.21) 


(2.22) 


The  corresponding  stress  increments  required  by  £qs.  (2.1), 
(2.2)  and  (2.3)  are 


>n+1  = pn  + * Y.  ^tn+J|5 


(2.23) 
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(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 


Equation  (2.17)  gives  the  total  volumetric  strain 
rate.  Equations (2.18)  through  (2.22)  are  the  total  deviatoric 
strain  rates.  The  isotropic  component  of  the  stress  tensor 
is  given  in  Eq.  (2.23).  Note  that  this  component  is  nega- 
tive in  compression,  i.e.,  it  is  opposite  in  sign  from  the 
hydrodynamic  pressure. 


2.6  NONLINEAR  MATERIAL  BEHAVIOR 

The  basic  assumption  of  nonlinear  continuum  mechanics 
is  that  the  total  strain  rate  e. . is  composed  of  elastic  n. . 

_ , t 1 j ^ *1 

and  inelastic  tt„  components,  i.e.,  J 

eij  = nij  + ^ij-  (2.29) 

Only  the  elastic  components  are  used  to  increment  stress, 
if  "jj  1*  0,  then  all  total  strain  rates  in  Eqs.  (2.23) 
through  (2.28)  must  be  replaced  by  elastic  strain  rates.  The 
constitutive  equation  determines  at  a given  stress  state. 

Various  models  for  nonlinear  material  behavior  have 
been  developed  and  used  at  S3  in  one  and  two  dimensional 
stress  wave  codes.  These  include  models  of  tension  failure, 
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plastic  flow,  effective  stress,  irreversible  pore  collapse 
and  dilatancy.  At  the  present  time,  TRES  contains  only  the 
plastic  flow  model,  which  has  been  used  to  simulate  non- 
linear material  behavior  in  the  fault  zone. 

For  plastic  flow  we  assume  that  the  maximum  allowable 

A 

stress  difference,  Y,  is  known.  The  stress  deviators 
are  calculated  using  total  strain  rates  as  given  by  Eqs. 
(2.24)  through  (2.28),  where 


S 

1 1 


9 


S 

22 


S 

1 3 


a 


xz 


t 


S 

2 3 


The  second  deviatoric  invariant  is 


J' 
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*1  A A 

1 s.  .s. . 

2 ij  ji 
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a + a a + a + c * + a2 

y x y xy  xz  yz 


Plastic  flow  results  if 


J' 

2 


(2.30) 


Then 


(2.31) 


and  the  stresses  are  used  to  calculate  accelerations. 
The  increment  in  plastic  strain  corresponding  to  the  stress 
adjustment  given  by  Eg.  (2.31)  is 


/s. 


(2.32) 
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The  corresponding  increment  in  plastic  work  is 


AE  = S . . Arr . . 
ID  ID 


(2.33) 


A plastic  work  rupture  criterion  has  proven  to  be  an  effec- 
tive technique  for  controlling  the  rupture  velocity  during 
earthquake  simulation  (Cherry,  et  aL,  1976). 

The  flexibility  of  the  code  to  accept  various  modes 
of  nonlinear  material  behavior  is  contained  in  the  stress 
rate-strain  rate  formulation  given  by  Eqs.  (2.23)  through 
(2.28).  Modification  of  these  stresses  at  time  tn+^  is 
obtained  through  the  constitutive  equation  which  determines 
the  inelastic  strain  rate  in  Eq.  (2.29).  The  code  now  con- 
tains the  plastic  flow  model.  Additional  constitutive 
models  are  a straightforward  extension  of  those  existing 
in  one  and  two  dimensional  codes.  These  will  be  included 
as  the  need  arises. 


2.7  THE  LINEAR  MODULE 

The  nonlinear  region  in  TRES  is  linked  to  a linear 
module.  The  physical  dimensions  of  the  surface  over  which 
the  nonlinear  and  linear  modules  are  linked  must  be  speci- 
fied by  the  user. 

The  link  has  two  objectives,  decreasing  computation 
time  and  coupling  the  results  of  the  finite  difference 
calculation  to  analytic  wave  propagation  techniques. 

The  former  is  achieved  since  stress  is  related  to 
s^rain  through  Hooke's  law.  Stress  is  not  a saved  variable 
in  the  linear  module  since  it  is  calculated  from  nodal 
point  displacements. 

The  latter  is  achieved  by  constructing  an  equivalent 
point  source  which  provides  an  analytic  continuation  of  the 
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outgoing  stress  field  in  the  elastic  region.  The  divergence 
and  curl  of  the  displacement  field  are  monitored  over  a 
spherical  surface  in  TRES  and  are  used  to  obtain  the 
multipole  coefficients  in  a spherical  wave  expansion  of  the 
radiating  elastic  field. 

In  the  following  section  we  describe  the  application 
of  TRES  for  simulating  strike-slip  earthquake  faulting. 

The  determination  of  an  equivalent  elastic  point  source 

• 

representation  of  the  earthquake  from  the  TRES  output  in 
the  linear  regime  is  an  important  step  in  describing  the 
character  of  the  source  as  a radiator  of  far-field  seismic 


waves. 


III.  FINITE  DIFFERENCE  EARTHQUAKE  MODEL 
3.1  INTRODUCTION 

It  is  only  within  the  last  10-15  years  that  serious 
attempts  have  been  made  to  develop  a mathematical  model  for 
the  earthquake  source.  Source  theory,  as  it  is  often  called, 
remains  a subject  of  considerable  interest  and  controversy 
and  much  scientific  energy  is  being. expended  to  develop  new 
source  theories,  expand  old  ones  or  to  demonstrate  the  abi- 
lity of  a particular  theory  to  explain  some  portion  of 
earthquake  data.  The  subject  of  this  section  of  the  report 
is  an  earthquake  source  model  that  employs  finite  difference 
numerical  techniques.  We  use  chis  model  to  carry  out  a 
three-dimensional  simulation  of  a small  earthquake. 

Most  of  the  source  models  discussed  in  the  literature 
are  elastodynamic  models.  That  is,  the  earthquake  source 
is  characterized  as  a radiator  of  elastic  waves  and  the  non- 
linear processes  near  the  fault  are  not  explicitly  treated. 
The  formulation  of  these  models  is  essentially  analytical, 
though  machine  calculations  may  be  required  to  evaluate 
some  of  the  expressions.  On  the  other  hand,  numerical 
techniques  (finite  difference  or  finite  element)  have  the 
potential  capability  to  model  fine  details  of  the  rupture 
process.  There  are  at  least  three  important  applications 
for  which  such  detailed  models  can  be  used : 

1.  As  a research  tool,  to  develop  an  understanding 
of  the  details  of  the  rupture  process. 

2.  As  a research  tool,  to  verify  and  extend  the 
assumptions  and  approximations  made  in  develop- 
ing analytical  elastodynamic  models. 

3.  To  numerically  simulate  earthquake  ground  motion 
in  the  near-field,  strong  motion  regime  where 
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analytical  approximations  are  unavailable  or 
unapplicable.  The  applications  are  either  for 
research  or  directly  to  earthquake  engineering 
problems  of  immediate  interest. 

There  have  been  a few  attempts  to  simulate  earthquake 
faulting  with  numerical  techniques.  The  work  of  Dieterich 
(1973) , Andrews  (1975) , and  Geller  and  Frazier  (1976)  is 
known  to  us.  Cherry,  et  al . , (1976)  incorporated  a two- 

dimensional  (plane  strain)  stick-slip  rupture  model  into  a 
finite  difference  code.  The  model  used  here  is  essentially 
a three-dimensional  extension  of  this  earlier  model. 

The  results  presented  in  this  report  may  be  viewed  in 
three  parts.  First,  we  describe  the  numerical  earthquake 
simulation,  showing  that  plausible  and  intuitively  appealing 
behavior  occurs  in  the  fault  vicinity.  Second,  the  earthquake 
generated  elastic  wave  field  is  decomposed  into  a unique  and 
exact  (in  a homogeneous  space)  equivalent  elastic  point 
source.  The  ability  to  obtain  such  an  equivalent  source 
from  the  results  of  a finite  difference  calculation  is  an 
important  step  because  it  allows  a linkage  between  the  numeri- 
cal codes  and  the  analytical  elastic  wave  propagation  techni- 
ques of  theoretical  seismology.  Finally,  we  compare  our 
earthquake  source  to  equivalent  events  from  volume-  and  dis- 
location-type elastodynamic  models.  The  interesting  result 
is  that  as  a radiator  of  elastic  waves  our  source  is  much 
more  like  the  volume  source  model. 

3.2  MODELING  PROCEDURE 

In  order  to  develop  a numerical  model  for  simulating 
the  earthquake  faulting  process  it  is  necessary  to  specify: 

1.  The  region  in  which  faulting  is  to  occur. 

The  constitutive  behavior  of  the  material  in  the 
vicinity  of  the  fault. 


2. 


3.  The  condition  for  rupture  initiation. 

4.  The  criterion  controlling  adjustment  of  the 
stress  field  during  rupture. 

5.  A criterion  allowing  the  rupture  to  heal. 

All  earthquake  models  must,  implicitly  or  explicitly,  deal 
with  these  five  items. 

Finite  difference  techniques  allow  almost  unlimited 
freedom  for  modeling  earthquake  faulting  as  outlined  above. 
The  restrictions  are  generally  economic  rather  than  techni- 
cal. In  this  section  we  wish  to  describe  a specific  model 
we  chose  to  study,  but  also  to  point  out  the  generality  of 
the  method  and  the  alternative  choices  that  can  be  made. 

A three-dimensional  finite  difference  earthquake 
model  was  developed  based  on  experience  with  two-dimensional 
modeling  studies.  The  basic  difference  equations  are  simi- 
lar to  those  reported  by  Cherry,  et  al^_  (1970).  The  stick- 
slip  rupture  model  is  an  extension  of  that  used  by  Cherry, 
al • (1976),  in  a study  of  the  ground  motion  from  two- 
dimensional  (plane  strain)  earthquake  faulting.  The  numeri- 
cal scheme  for  solving  the  differenced  form  of  the  governing 
differential  equations  incorporates  explicit  time  stepping 
for  updating  the  stress.  Details  have  been  presented  in 
Section  II. 

The  geometry  for  the  calculation  is  indicated  in 
Figure  3.1.  A rectangular  fault,  1.0  x 0.6  km2,  was  speci- 
fied as  the  surface  on  which  slip  was  to  be  allowed.  That 
is,  a discontinuity  in  X-displacement  was  allowed  across 
this  surface.  It  is  clear  that  some  nonlinear  mechanism  must 
be  postulated  for  the  region  adjacent  to  the  slip  surface. 

In  fact,  on  the  edges  of  this  surface  the  strain  is  unbounded. 
The  questions  that  must  be  addressed  are:  What  is  the  non- 

linear material  behavior  in  the  fault  zone?  What  is  the 
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extent  of  the  fault  zone;  that  is,  how  large  is  the  region 
in  which  material  nonlinearities  are  important? 

The  fault  zone  was  specified  to  be  a rectangular 
parallelepiped  2.0  x i.6  x 0.8  km3  as  indicated  in  Figure 
3.1.  Outside  this  zone  the  material  behavior  was  linear 
elastic.  The  elastic  P and  S wave  velocities  are  a = 5.7 
km/sec  and  0 = 3.4  km/sec.  The  density  is  p = 2.8  gm/cm3 
and  the  shear  modulus  is  p = 324  kbar.  In  the  fault  zone, 
elastic-plastic  material  response  consistent  with  a von 
Mises  yield  criterion  was  allowed.  That  is,  if  the  second 

deviatoric  invariant  (J')  is  greater  than  a specified  value, 
then 


2 


A 

where  is  the  stress  deviator  tensor,  S.^ . is  the  adjusted 
stress  deviator  and  . 13 


J' 

? 


(3.2) 


For  a triaxial  test  Y corresponds  to  the  maximum  stress 
difference  at  failure. 

Another  important  feature  of  the  calculation  are  the 
grid  dimensions  for  the  finite  difference  mesh.  All  zones 
were  cubes,  100  meters  on  a side.  Then  the  fault  surface 
itself  is  6 x 10  grid  dimensions  while  the  fault  zone  in- 
cludes 20  x 16  x 8 grid  dimensions.  With  a grid  this  coarse 
stress  waves  with  frequencies  greater  than  3-5  Hz  are  not 
propagated  accurately  and  this  must  be  considered  when 
evaluating  the  results.  In  order  to  obtain  an  equivalent 
elastic  source  (described  in  a subsequent  section) , it  was 
necessary  to  propagate  the  full  time  history  of  the  earth- 
quake generated  ground  motion  to  some  elastic  radius.  This 
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was  selected  to  be  at  1.5  km  as  indicated  in  Figure  3.1. 

The  time  step  (which  is  controlled  by  the  grid  dimensions) 
was  0.005  seconds  and  225  time  steps  were  required  to  com- 
plete the  problem.  Taking  advantage  of  the  symmetry  of  the 
problem,  the  calculation  was  actually  run  for  only  a quadrant 
of  the  whole  space,  that  where  X and  Z are  positive.  A de- 
tailed description  of  the  symmetry  conditions  imposed  on  the 
boundary  of  the  finite  difference  grid  appears  in  Savino, 
et  al.  (1975) , Appendix  A. 

The  elastic-plastic  stress-strain  behavior  is  one 
of  the  simplest  nonlinear  mechanisms  that  can  be  introduced 
to  dissipate  the  stress  concentrations  at  the  edges  of  the 
dislocation.  Certainly,  the  size  of  the  fault  zone  and  the 
constitutive  behavior  of  the  material  within  this  zone  are 
important  parameters  of  the  calculation.  While  the  finite 
difference  formulation  allows  great  freedom  in  how  they  are 
specified,  we  made  this  relatively  simple  choice  as  being 
most  appropriate  for  our  first  three-dimensional  calculation. 

How  is  rupture  to  be  initiated?  In  the  two-dimen- 
sional earthquake  simulations  by  Cherry,  et  al.  (1976), 
rupture  initiation  was  controlled  by  a criterion  based  on 
the  amount  of  plastic  work  accumulated  at  a location  on  the 
fault  surface.  However,  the  rupture  velocity  was  found  co 
be  nearly  constant.  Then  a simpler  though  nearly  equivalent 
crite^ion  is  to  specify  the  rupture  velocity  and  this  was 
the  procedure  followed  in  the  calculation  described  here. 

The  rupture  was  assumed  to  grow  symmetrically  with  velocity 
Vr  = 2.7  km/sec  (80  percent  of  the  shear  wave  velocity)  un- 
til reaching  the  boundary  of  the  specified  fault  surface. 

The  fault  growth  pattern  is  indicated  in  Figure  3.1.  The 

prestress  was  assumed  to  be  pure  shear:  cr^  =1.0  kbar. 

— 

Therefore,  in  the  coordinate  system  of  Figure  3.1,  the 
faulting  is  bilateral  and  the  event  is  right-lateral  strike- 
slip. 
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The  next  question  to  be  answered  is:  How  does  the 

stress  field  adjust  during  rupture?  Initially,  all  the 
material  in  the  fault  zone  was  specified  to  be  on  the  fault 
surface.  In  view  ^f  (3.1)  this  means  Y = /3J>  or  Y = /3 
kbar  since  at  the  beginning  of  the  calculation  J'  . a(°)2  = 

1.0  kbar2.  Then  all  the  material  in  the  fault  zone  flows 
plastically  when  J'  increases  beyond  its  initial  value.  How- 
ever, this  increase  in  J',  leading  to  plastic  flow,  will  be 
restricted  to  the  region  around  the  tip  of  the  fault  since 
elsewhere  the  effect  of  the  rupture  is  to  relax  the  shear 
stress  from  its  initial  value.  At  the  time  of  rupture 
initiation,  grid  points  on  opposite  sides  of  the  fault  sur- 
face are  uncoupled  and  allowed  to  slip  as  the  tangential 
stress  is  relaxed  (over  three  time  steps)  to  its  kinetic 
friction  value  which  was  taken  to  be  0.85  kbar.  The  nominal 
dynamic  stress  drop  is  then  Aa  = 1.0  - 0.85  kbar  - 150  bars. 

A given  location  on  the  fault  sticks  when  there  is  a rever- 
sal in  relative  velocity  Between  points  on  opposite  sides 
of  the  fault. 

The  tangential  stress  at  several  locations  on  the 
fault  surface  is  plotted  in  Figure  3.2  at  nine  time  points 
at  which  the  grid  was  monitored.  We  see  that  the  stress 
can  build  up  slightly  before  rupture  initiation  at  which 
time  it  is  relaxed  to  0.85  kbar.  After  the  fault  is  tied 
the  stress  continues  to  change  before  reaching  its  static 
value.  The  static  values  vary  over  the  fault  surface  from 
0.77  - 0.91  kbar  with  the  average  value  being  0.854  kbar. 
Therefore,  the  average  static  stress  drop  is  146  bars, 
quite  close  to  the  nominal  Act  of  150  bars. 

In  Figure  3.3  is  plotted  the  X-displacement  corresponding 
to  the  tangential  stress  of  Figure  3.2.  Once  again,  the  sampl- 
ing is  only  at  nine  evenly  spaced  time  points.  However,  we 
do  know  the  time  of  rupture  initiation  at  each  point.  Taking 
this  into  account,  we  see  that  the  displacement  time  history 
on  the  fault  is  consistent  with  a simple  ramp  with  a rise 
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time  of  0.08  - 0.10  seconds,  depending  on  position.  This 
interpretation  has  been  used  to  draw  the  curves  in  the 
figure.  We  caution  that  we  do  not  mean  to  infer  that  these 
curves  show  the  actual  time  history  of  displacement,  only 
that  the  known  values  are  consistent  with  this  simple  model. 

beginning  of  this  section  we  listed  five  items 
that  must  be  specified  for  a deterministic  earthquake  fault- 
ing model.  In  specifying  these  items  we  have  chosen  to  keep 
the  problem  as  simple  as  possible,  consistent  with  the 
physics  as  we  understand  it.  The  result  is  a rupture  pro- 
cess characterized  by  a rather  complex  adjustment  of  the 
tangential  stress  across  the  fault  (Figure  3.2),  but  a 
relatively  simple  time  history  of  displacement  on  the  fault 
surface.  In  the  way  "dislocation"  models  of  earthquake 
faulting  are  usually  formulated  (e.g.,  Haskell,  1964),  dis- 
placement time  histories  very  much  like  those  in  Figure  3.3 
are  specified  on  the  fault  surface.  In  using  such  models, 

attention  is  generally  given  to  the  time  history  of 
shear  stress  on  the  fault.  In  our  model  the  adjustment  of 
the  shear  stress  (Figure  3.2)  is  clearly  influenced  by  the 
material  behavior  m the  fault  zone. 

As  pointed  out  above,  the  rupture  process  in  our 
finite  difference  model  is  rather  similar  to  that  assumed 
in  constructing  the  dislocation-type  models  of  earthquake 
faulting  of  which  there  are  many  examples  in  the  literature 
(e.g.,  Haskell,  1964;  Savage,  1966;  Brune,  1970;  Dahlen, 

1974;  Richards,  1976).  One  difference  is  the  plastic  flow 
allowed  in  the  fault  zone  to  relieve  the  stress  concentra- 
tions that  occur  at  the  crack  edges.  Another  difference  is 
that  the  fault  surface  does  not  remain  plane,  but  is  warped 
as  is  seen  in  Figure  3.4.  There  the  out-of-plane  (Y)  dis- 
placement of  the  fault  is  shown  versus  position  at  tl?e  last 
time  step  and  at  nine  snapshots  in  time  for  the  displace- 
ments on  the  X-axis.  The  maximum  out-of-plane  displace- 
ments are  at  the  end  of  the  fault  where  they  are  as  much  as 
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5 cm  or  approximately  20  percent  of  the  average  fault  off- 
set. 

Our  primary  interest  in  this  paper  is  in  the  source 
as  a radiator  of  seismic  waves.  As  a first  indication  of 
the  character  of  the  radiated  field,  some  typical  velocity 
and  displacement  time  histories  are  shown  in  Figure  3.5  at 
two  stations  1 . 5 km  from  the  center  of  the  fault.  In  these 
time  histories  much  of  the  complexity  can  be  attributed  to 
the  fact  that  near-  and  far-field  P and  S wave  effects  are 
all  represented.  In  the  next  section  the  outgoing  wave 
field  is  represented  in  equivalent  elastic  source  form  from 
which  it  can  be  decomposed  into  itrj  various  components. 


3.3  EQUIVALENT  ELASTIC  SOURCE  REPRESENTATION 

Finite  difference  methods  are  useful  for  modeling 
details  of  the  rupture  physics  but  are  not  practical  for 
propagating  the  resulting  seismic  waves  very  far  from  the 
-fault.  In  the  seismic  wave  (linear  elastic  wave  propagation) 
regime  the  analytical  techniques  of  theoretical  seismology 
are  generally  much  more  efficient  and  economical.  In  order 
to  bridge  the  gap  between  complex  source  calculations  such 
as  that  described  in  the  previous  section  and  elastic  wave 
propagation  theories,  it  is  necessary  to  construct  an  equi- 
valent elastic  source.  As  was  pointed  out  by  Bache  and 
Harkrider  (1976) , a convenient  equivalent  source  can  be  ob- 
tained by  expanding  the  outgoing  elastic  wave  field  in 
spherical  harmonics.  The  procedure  is  briefly  outlined  be- 
low. 

The  Fourier  transformed  equations  of  motion  in  a homo- 
geneous, isotropic,  linear  elastic  medium  may  be  written 
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where  u is  displacement  and  ka  and  kg  are  the  congressional 
and  shear  wave  numbers.  The  Cartesian  potentials  x(4>  and 
X are  defined  by 


_ , (3.4) 

X 2 7 x S • 

The  potentials  can  be  shown  to  satisfy  wave  equations  and 

therefore  can  be  expanded  in  spherical  eigenfunctions  as 
follows : 


00  $, 

X ^ (R,w)  = h|2^  (k jR)  |a|^  ^ cos 

1=0  m=0  (3.5) 

+ B^’  (hi)  sin  m*l  p”(cos8),  j = 1,2, 3, 4 

k,  " ka  " and  ki  = ks  = m/8  for  i > 1,2,3.  The 
h£  are  spherical  Hankel  functions  of  the  second  kind  and 
the  are  associated  Legendre  functions.  The  vector  R 
has  as  components  the  usual  spherical  coordinates  R,0,$, 
referred  to  the  Cartesian  system  of  Figure  3.1. 

Equations  (3.5),  together  with  (3.3),  provide  an 
elastic  point  source  representation  of  the  (outgoing)  dis- 
placement field.  The  values  of  the  multipole  coefficients, 

Alm  (lo)'  B*m  ' 3 = If 2, 3, 4,  prescribe  the  displacement 
field  at  all  points  in  the  homogeneous  medium  where  (3.3) 
applies.  Therefore,  this  is  a convenient  source  form  for 
studying  whole  space  radiation  properties,  decomposing  the 
field  into  near-  and  far-field  components,  etc.  Minster 
(1973)  gives  an  extensive  discussion  of  the  manipulation  of 
sources  given  in  this  form  including  formulas  for  rotating 
the  source  to  any  desired  orientation  with  respect  to  a 
fixed  coordinate  system.  Further,  Harkrider  and  Archambeau 


(1976)  and  Bache  and  Harkrider  (1976)  have  derived  the  theory 
for  computing  surface  and  body  waves  in  layered  media  for 
sources  given  in  terms  of  multipole  coefficients. 

Using  the  orthogonality  of  the  spherical  harmonics, 
the  multipole  coefficients  are  related  to  the  potentials, 
XU\  by 
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and  R is  a radius  defining  a spherical  surface  on  which  the 
integral  (3.6)  is  to  be  computed.  For  the  finite  difference 
earthquake  calculation  discussed  here,  the  divergence  and 
curl  of  the  displacement  field  were  monitoried  at  R » 1.5  km 
the  "elastic  radius"  indicated  in  Figure  3.1.  since  a 
Cartesian  grid  was  used,  the  actual  implementation  required 
a linear  interpolation  to  the  spherical  surface  from  nearby 
grid  points.  Details  of  the  algorithm  for  computing  the 
integrals  (3.6)  from  finite  difference  code  output  are  given 
in  Appendix  G of  Bache,  et  al^  (1975).  The  tractability  of 
the  surface  integrand  in  (3.6)  is  clearly  dependent  on  the 
behavior  of  the  curl  and  divergence  as  functions  of  0 and 

4>.  For  this  earthquake  calculation  these  quantities  were 
quite  smooth  and  well-behaved. 
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Due  to  the  symmetry  of  the  source,  only  even  order 
U = 0,2,3,...)  multipoles  are  nonzero.  The  odd  order 
terms  are  associated  with  unilateral  propagation  effects. 
The  1=2  terms  represent  the  quadrupole  or  double-couple 
radiation  field  generated  by  the  source  while  the  higher 
order  terms  are  perturbations  on  the  double-couple  field 
introduced  by  propagation  effects.  Also,  there  is  a single 
monopole  term,  (w) , which  represents  a spherically  sym- 
metric rotation  about  the  Z-axis.  The  monopole  shear  term 
arises  as  a consequence  of  the  out-of-plane  rotation  of  the 
fault  surface  (Figure  3.4).  As  expected  from  the  require- 
ment for  the  source  region  to  reach  static  equilibrium,  this 
term  vanishes  at  long  time. 

The  amplitude  of  the  monooole  and  largest  of  the 
quadrupole  and  octupole  coefficients  is  plotted  in  Figure 
3.6.  Note  that  at  low  frequencies  the  source  is  dominated 
by  the  quadrupole  terms.  In  fact  these  terms  fall  off  with 
w2  at  low  frequencies  while  terms  of  other  orders  decay 
with  higher  powers  of  w.  At  higher  frequencies  more  and 
more  terms  are  required  to  insure  convergence  of  (3.5). 

3 *3  4 THE  seismic  wave  radiation  character  of  the  source 

Having  the  multipole  coefficients,  computation  of 
displacement  spectra  at  any  location,  R,  is  a straightfor- 
ward application  of  (3.3)  and  (3.5).  Of  course,  it  is  neces- 
sary to  compute  enough  terms  so  that  (3.5)  converges.  In 
the  sequel  we  will  restrict  attention  to  the  far-field  por- 
tion of  the  radiation  field  which  is  defined  as  that  portion 
which  decays  as  R"1.  Then  the  "P«  and  "S"  wave  portions  of 
the  wave  field  are  conveniently  separated  with  the  first 
term  of  (3.3)  giving  the  "P  wave”  and  the  second  term  giving 
the  "S  wave".  In  fact,  it  can  easily  be  shown  that: 


47 


Frequency  (Hz) 


Figure  3.6.  The  amplitude  of  the  monopole  and  largest  quadru— 
pole  and  octupole  terms  for  the  equivalent  elastic 
source  representation  of  our  earthquake  referred 
to  the  coordinate  system  of  Figure  3.1. 
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and  the  derivative  with  respect  to  R reduces  to 
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The  geographic  coordinate  system  and  fault  gee 
for  discussing  the  far-field  stress  wave  radic--  ^n  pro,__.r- 
ties  of  the  source  are  shown  in  Figure  3.7.  * -sition  with 

respect  to  the  fault  is  denoted  by  a ray  (A)  specified  by 
its  azimuth  (<f>)  and  takeoff  angle  (t)  , where  t = it  - 0 . 

The  radial  distance  along  A is  denoted  by  R.  With  the  co- 
ordinate system  oriented  as  in  Figure  3.7f  the  u , u.  and 
__  R U 

u^  of  (3.7)  represent  the  P,  SV  and  SH  waves,  respectively. 

Far-field  displacement  spectra  are  shown  in  Figure 
3.8  at  two  positions  on  the  radiation  pattern.  The  spectra 
were  computed  at  R = 1.5  km  but,  in  view  of  (3.8),  can  be 
scaled  to  any  desired  range.  Terms  through  l = 4 were  re- 
tained in  the  multipolar  expansion.  These  are  sufficient  to 
insure  convergence  for  the  frequencies  plotted.  Note  that 
the  <j)  = 45° , t = 45°  position  is  the  same  as  the  left  plot 
in  Figure  3.5  where  total  radial  displacement  and  velocity 
were  plotted. 
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Figure  3.7.  Orientation  of  the  fault  with  respect  to  the  geo 
graphic  coordinate  system. 
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Figure  3.8.  Far-field  displacement  spectra  at  two  locations  (specified  by  4>  and  x)  on 
the  elastic  radius  (R  = 1.5  km). 


The  spectra  shown  in  Figure  3.8  are  remarkably  simple 
and  can  be  characterized  by  three  parameters;  the  long  period 
level,  the  corner  frequency  and  the  logarithmic  slope  of  the 
high  frequency  values.  The  long  period  level  fixes  the  seis- 
mic moment  for  the  event,  a quantity  that  will  be  discussed 
in  a later  section.  For  the  high  frequency  portion  of  the 
spectrum,  recall  that  the  grid  dimensions  impose  the  restric- 
tion that  frequencies  greater  than  3-5  Hz  are  not*  propagated 
accurately.  The  highest  frequencies  are  then  contaminated 
to  some  degree  by  these  effects  which  generally  tend  to 
increase  the  slope.  Also,  for  these  frequencies  the  spectra 
may  still  be  in  a transition  region  and  one  should  take  care 
when  comparing  to  true  high  frequency  asymptotes  obtained 
from  observations  or  analytic  theories. 

Corner  frequency  is  a somewhat  ambiguous  quantity 
with  a value  that  depends  on  the  measuring  technique.  Blandford 
(1975)  has  suggested  that  the  corner  frequency  be  defined  as 
that  freqiency  at  which  the  spectral  value  is  half  its  low 
frequency  limit.  This  definition  allows  a relatively  pre- 
cise determination  of  corner  frequency  and  will  be  used  here. 
Corner  frequency  is  radiation  pattern  dependent  as  is  indi- 
cated by  the  values  given  in  Table  3.1.  Note  that  for  fixed 
x the  corner  frequency,  P or  S,  decreases  as  the  observer 
moves  from  a position  normal  to  the  fault  ($  = 90°)  toward 
the  direction  of  rupture  propagation  (<J>  = 0°).  The  corner 
frequencies  of  Table  3.1  are  all  from  one  octant  of  the  focal 
sphere,  a sufficient  sampling  because  of  the  source  symmetry 
Taking  the  average,  the  corner  frequencies  are  1.38  for  the 

P wave  and  1.20  for  the  S wave  (SV  or  SH)  with  the  deviations 

from  these  values  being  less  than  0.1  Hz  in  all  but  a few 

cases.  On  the  average,  the  P wave  corner  frequency  is  1.15 

times  the  S wave  corner  frequency. 

The  corner  frequency  is  clearly  associated  in  some 
way  with  the  fault  dimensions.  Knowledge  of  the  relationship 
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between  the  two  affords  an  opportunity  to  obtain  estimates 
of  the  fault  dimensions  from  relatively  sparse  amounts  of 
data.  The  usual  assumption  is  that 


/ 

P/S 


(3.9) 


where  D is  the  fault  dimension,  f is  the  corner  frequency 
and  Cp,s  1S  a constant  for  P or  S waves.  For  example,  Brune 
(1970)  suggests  Cg  = 0.37.  For  our  computed  source  we  find 
from  the  values  of  Table  3.1,  using  D = 1.0  km,  that 

Cp  = 0.39  - 0.42, 

(3.10) 

Cg  * 0.30  - 0.46. 


Note  that  the  Cg  = 0.46  is  an  isolated  point  (SH  waves  at 
0 = t = 90°)  and  the  second  highest  value  is  C = 0.39.  The 

average  values  are  C = 0.41  and  C =0.35. 

P s 

From  Figure  3.8  it  seems  clear  that,  no  matter  how 
they  are  determined,  the  average  ratio  of  P to  S corner 
frequencies  (fp/fg)  is  greater  than  1 for  our  model.  This 
is  consistent  with  the  results  from  volume  source  models 
of  earthquakes  (e.g. , Archambeau,  1968).  We  saw  in  the  pre- 
vious section,  especially  Figure  3.3,  that  in  many  respects 
the  rupture  process  studied  here  behaves  very  much  like  that 
m a dislocation  model.  Peppin  and  Simila  (1976)  summarize 
the  corner  frequency  characteristics  of  many  of  the  proposed 
dislocation  models.  These  authors  point  out  that  disloca- 
tion models  with  subsonic  and  uniform  rupture  propagation 
(like  in  our  model)  give  fp/fg  < 1 (Savage,  1972;  Dahlen, 
1974).  On  the  other  hand,  those  dislocation  models  that 
have  fp/fs  > 1 either  have  supersonic  rupture  velocity 
(Brune,  1970;  Molnar , et  al. , 1973;  Burridge,  1975),  variable 
rupture  velocity  (Berckhemer  and  Jacob,  1968)  or  slip  con- 
centrated near  the  center  of  the  fault  (Sato  and  Hirasawa, 
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1973;  Madariaga,  1976).  Why  does  our  model  exhibit  corner 
frequency  behavior  more  like  these  latter  dislocation  models 
and  volume  source  models  than  a dislocation  model  such  as 
that  of  Savage  which  is  superficially  more  similar?  We  will 
return  to  this  question  in  later  sections  where  more  exten- 
sive comparison  to  analytical  models  is  made. 

We  complete  our  description  of  the  far-field  dis- 
placements from  our  finite  difference  earthquake  source 
with  the  radiation  patterns  shown  in  Figure  3.9.  Patterns 
are  computed  from  (3.7)  for  three  frequencies  and  two  levels 
of  truncation  of  the  sum  (3.5).  In  the  top  row  only  the 
monopole  and  quadrupole  terms  were  retained.  The  octupole 
term  was  added  for  the  patterns  of  the  bottom  row.  Only 
the  P and  SH  components  are  nonzero  at  6 = 90°.  We  observe 
that  at  low  frequency  the  source  is  nearly  pure  double 
couple,  with  the  only  variation  at  0.5  Hz  being  that  slightly 
more  SH  waves  are  propagated  normal  to  the  fault  (along  YQ) 
than  in  the  fault  direction  (along  XG) . At  higher  frequencie 
the  monopole  and  higher  order  terms  have  greater  effect, 
distorting  both  the  P and  SH  patterns.  Even  for  this  small 
and  rather  simple  event,  terms  other  than  the  quadrupole 
have  a significant  effect  on  the  short  period  (<^  1 sec) 
energy  radiated  to  the  far-field. 

The  propagation  of  substantially  larger  SH  waves  to 
azimuths  normal  to  the  fault  compared  to  azimuths  along  the 
fault  (e.g.,  50  percent  larger  at  1 Hz)  is  an  interesting 
result.  For  the  frequencies  plotted,  both  the  monopole  and 
octupole  act  in  this  direction.  To  the  extent  that  the 
rronopole  is  responsible,  the  out-of-plane  displacement  of  the 
fault  surface  is  a major  cause.  However,  the  same  thing 
can  occur  in  models  which  have  no  such  fault  surface  warping 
as  will  be  seen  in  the  following  section  for  an  analytical 
model.  In  this  case  the  cause  is  interference  between  waves 
from  either  end  of  the  bilateral  fault. 
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COMPARISON  TO  THE  ARCHAiiBEAU /MINSTER  RELAXATION  SOURCE 
MODEL 

We  saw  in  a previous  section,  especially  Figure  3.3, 
that  in  terms  of  a kinematic  view  of  the  faulting,  our 
finite  difference  model  seems  similar  to  dislocation  models. 
On  the  other  hand,  the  faulting  process  is  initiated  by  the 
introduction  of  a weak  zone  in  a prestressed  medium  and  in 
that  sense  is  similar  to  a relaxation  model.  Further,  there 
is  a finite  volume  within  which  nonlinear  material  response 
is  allowed.  In  short,  our  finite  difference  source  bears 
superficial  similarities  to  both  relaxation-volume  source 
models  and  dislocation  models.  In  the  remainder  of  the 
paper  we  shall  compare  our  results  to  those  obtained  from 
various  forms  of  these  analytical  models. 

The  Archambeau/Minster  elastodynair.ic  source  theory 
(Archambeau,  1968;  Minster,  1973)  describes  the  earthquake 
in  terms  of  the  relaxation  of  a pre-existing  stress  field 
about  a zone  of  weakness  in  which  the  rigidity  modulus 
vanishes.  The  theory  has  been  cast  in  a computationally 
convenient  form  for  several  geometries  in  which  the  weak 
zone  is  formed  by  a growing  and/or  propagating  spherical 
volume.  The  displacement  field  solutions  for  the  Archambeau/ 
Minster  source  theory  are  given  directly  in  terms  of  the 
multipolar  expansion  (3.5).  Since  we  have  described  our 
finite  difference  earthquake  source  in  the  same  form,  com- 
parison between  the  two  is  conveniently  accomplished. 

Formulation  of  the  Archambeau/Minster  source  model 
requires  specification  of  o and  0,  the  elastic  wave  veloci- 
ties; VR,  the  rupture  velocity;  and  S(0*,  the  prestrain 
which  is  equal  to  the  ratio  of  prestress  (a ^ ) to  shear 
modulus  (y) . To  control  the  growth  of  the  weak  zone  we 
specify  two  lengths,  R(t)  and  d(t),  which  are  functions 
of  time.  The  R(t)  is  the  radius  of  the  spherical  volume 
of  vanishing  rigidity,  while  d (t)  is  the  distance  of  the 
center  of  this  sphere  from  the  origin.  Then  R(t)  + d (t) 


r" 


grows  steadily  with  velocity  VR.  The  parameter  R , intro- 
auced  by  Archambeau  to  represent  the  dimension  o^the  pre- 
stressed region,  will  oe  infinite  and  play  n c part  in  the 
solution. 

Three  source  geometries  are  of  interest: 

1.  Stationary  expanding  sphere:  if  we  set  d(t)  “0, 

the  spherical  region  expands  uniformly  about  its 
center  point.  This  form  of  Archambeau' s model 
has  been  used  (Archambeau,  1972;  Archambeau  and 
Sammis,  1970;  Bache,  1976)  in  a number  of  studies 
of  the  tectonic  strain  release  associated  with 
large  underground  explosions. 

2.  Tangentially  expanding  sphere:  If  d (t)  = R(t), 

the  spherical  region  expands  from  a fixed  point 
on  the  boundary.  The  radiation  field  from  this 
geometry  is  extensively  discussed  by  Minster 
(1973). 

3.  Growing  and  propagating  sphere:  The  d(t)  and 

R(t)  can  take  any  functional  form  subject  to 

the  restriction  that  d(t)  + R(t)  = V t. 

R 

For  comparing  with  our  finite  difference  earthquake 
the  selection  of  appropriate  values  for  most  of  the  para- 
meters is  obvious.  We  take  a = 5.7  km/sec,  0 = 3.4  km/sec 
and  VR  = 2.7  km/sec.  Since  the  stress  drop  is  total  in  the 
analytic  model,  we  take  a(0)  = 150  bars  or  S(0)  = 4.63  * 

10  4.  Specification  of  the  weak  zone  dimensions  depends  on 
the  model.  Also,  we  must  account  for  the  fact  that  the 
finite  difference  fault  is  bilateral. 

The  simplest  comparison  is  made  by  using  the  stationary 
expanding  sphere,  model  1 above,  and  associating  the  diameter 
with  the  fault  length  (L)  of  1.0  km.  That  is,  let  R(t) 
grow  with  velocity  VR  to  a maximum  of  0.5  km.  This  source 
approximates  our  bilateral  fault  in  the  sense  that  it  is 
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growing  equally  in  all  directions.  However,  the  source  is 
a pure  double-couple;  there  are  no  higher  order  terms  repre- 
senting the  interference  between  waves  originating  in  dif- 
ferent regions  of  the  fault. 

In  Figure  3.8  were  shown  far-field  displacement  spectra 
for  our  finite  difference  earthquake.  The  same  quantities 
were  computed  for  the  stationary  expanding  sphere  model 
described  above  and  the  two  are  compared  for  x = 30°, 

$ = 150 > in  the  left  plot  in  Figure  3.10  which  is  typical 
of  the  entire  radiation  pattern. 

The  results  of  this  first  comparison  are  simply  sum- 
marized: The  long  period  level  and  relative  excitation  of 

P,  SV  and  SH  waves  agree  almost  exactly  while  the  corner 
frequency  is  considerably  higher  for  the  analytical  model. 

As  pointed  out  by  Minster  (1973),  the  long  period  level  for 

either  of  the  expanding  sphere  models  is  controlled  by  the 

(0)  2 
strain  drop  (S  ) and  the  volume  of  the  final  sphere  (=  L3). 

The  corner  frequency  depends  on  V_/L. 

K 

The  natural  question  to  ask  is:  How  much  is  the 

short  period  portion  of  the  spectrum  influenced  by  propaga- 
tion (interference)  effects?  The  tangentially  expanding 
sphere  is  the  simplest  Archambeau/Minster  model  that  includes 
the  effects  of  asymmetric  rupture  propagation.  To  simulate 
a bilateral  rupture  we  need  to  use  two  such  models,  one  with 
the  rupture  traveling  to  the  right  (SR)  and  one  with  it 
traveling  to  the  left  (S^) . To  maintain  the  long  period 
spectrum  at  the  proper  level  we  take  L = 1.0  km  in  both 
models  and  obtain  the  bilateral  fault  (S_)  from  S = (S  + 

SI /'*•  As  a result,  consists  of  the  even  order  multi- 
poles  (A  = 2,4,6,...)  from  SR  or  while  the  odd  order 
(A  = 3,5,7,...)  multipoles  drop  out. 

The  comparison  of  the  finite  difference  earthquake 
spectra  and  the  bilateral  model  described  above  is  shown 
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on  the  right  in  Figure  3.10.  The  agreement  between  long 
period  spectra  is  again  nearly  exact.  We  see  that  addi- 
tion of  the  higher  order  interference  terms  does  consider- 
ably enhance  the  agreement  in  the  short  period  region.  Sub- 
stantial disagreement  between  the  two  models  seems  to  be 
restricted  to  the  vicinity  of  the  corner  frequency  where 
the  analytical  model  breaks  much  more  sharply  to  its  high 
frequency  slope. 

Since  both  the  analytic  and  finite  difference  source 
models  are  given  in  multipolar  form,  the  most  detailed  com- 
parison is  via  the  multipole  coefficients.  To  the  extent 
that  these  agree,  the  entire  radiation  pattern,  near-  and 
far-field,  is  the  same.  Comparison  of  the  amplitudes  of 
the  multipoles  (for  the  fault  orientation  of  Figure  3.7) 
is  shown  in  Table  3.2  at  three  points  in  the  low,  intermediate 
and  high  frequency  range.  For  both  the  finite  difference 
model  (FD)  and  the  analytical  model  (TES)  we  show  the  values 
of  the  nonzero  coefficients,  through  the  octupole,  normalized 
to  the  leading  term,  A^  (FD)  , for  the  finite  difference 
model . 

From  Table  3.2  it  is  clear  that  at  low  frequency  the 
FD  source  is  essentially  a pure  double-couple  and  is  nearly 
identical  to  the  TES  source.  Agreement  between  the  two  is 
increasingly  worse  as  frequency  increases.  A notable  fea- 
ture of  the  high  frequency  FD  values  is  the  importance  of 

the  out-of-plane  fault  displacement  terms  A^  and  A^ 

0 0 2 0 

which  do  not  appear  in  the  analytical  solution.  From  this 
comparison  we  would  expect  the  radiation  patterns  for  the 
two  sources  to  be  substantially  the  same  at  low  frequencies 
and  differ  as  frequency  increases.  This  is  illustrated  in 
Figure  3.11  where  patterns  for  the  TES  that  are  directly 
comparable  to  the  FD  patterns  of  Figure  3.9  are  shown. 

The  availability  of  an  analytical  model  that  has  an 
equivalent  source  with  so  many  features  in  common  with  our 
FD  equivalent  source  is  quite  useful  in  attempting  to 
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Table  3.2 


. Comparison  of  Multipole  Coefficients  for  the 
Finite  Difference  Earthquake  (FD)  to  Those 
from  the  Analytical  Bilateral  Fault  Constructed 
from  the  Archambeau/Minster  Tangentially  Ex- 
panding Sphere  (TES)* 


Frequency  (Hz) 

0 

.1 

1 

.0 

5 

.0 

Aim’/An’  tFD> 

FD 

TES 

FD 

TES 

FD 

TES 

Monopole 

a(3) 

• 

0 0 

Quadrupole 

a'1’ 

2 1 

0.01 

0 

0.21 

0 

0.63 

0 

1.00 

0.97 

1.00 

0.89 

1.00 

0.42 

b'2> 
2 1 

1.01 

0.97 

1.13 

0.89 

0.67 

0.42 

a‘3> 
2 0 

0.06 

0 

0.32 

0 

0.84 

0 

a<3> 

22 

1.04 

0.97 

1.09 

0.89 

0.46 

0.42 

b<4> 

22 

0.16 

0.15 

0.19 

0.15 

0.19 

0.05 

Octupole 

A(1> 

4 ! 

- 

m 

0.02 

0.05 

0.27 

0.26 

A(1> 

4 3 

- 

- 

- 

- 

0.01 

0.02 

b(2> 

4 1 

- 

- 

0.03 

0.02 

0.27 

0.09 

b'2> 

4 3 

- 

- 

- 

- 

0.01 

0.02 

a<3> 

4 0 

- 

- 

0.11 

0.07 

0.75 

0.36 

a(3> 

42 

- 

- 

- 

0.01 

0.07 

0.06 

a<3> 

4 4 

- 

- 

- 

- 

0.01 

0.02 

b'4> 

42 

- 

- 

- 

_ i 

0.01 

0.02 

b(4) 

4 4 

* 

The  dashed  entries  signify  a nonzero  value  less  than  0.01. 
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Figure  3.11.  Far-field  radiation  patterns  for  the  bilateral  fault  (TES)  constru 
from  the  tangentially  expanding  sphere  model.  These  patterns  are 
6 = 90°  and  are  directly  comparable  to  those  in  Figure  3.9  for  the 
finite  difference  source.  Again,  a = 9/R  cm/km. 


understand  our  results,  isolate  the  origin  of  the  observed 
effects,  etc.  Therefore,  we  might  be  interested  in  perturb- 
ing the  Archambeau/Minster  model  to  get  an  equivalent  source 
that  is  even  closer  to  the  FD.  The  third  geometry  listed 
above,  the  growing  and  propagating  sphere  model,  allows  con- 
siderable freedom  to  adjust  the  radiated  wave  field  as  re- 
quired, though  this  will  not  be  pursued  further  here. 

To  summarize  the  results  of  this  comparison  between 
analytical  and  finite  difference  sources,  the  first  point  is 
that  when  the  fault  lengths  are  the  same,  the  two  sources 
are  nearly  identical  at  long  periods.  This  agreement  in 
seismic  moment  is  a remarkable  result  to  which  we  will  re- 
turn. Agreement  at  high  frequencies  is  less  impressive,  but 
much  can  be  learned  from  detailed  comparison  of  the  two 
equivalent  sources.  Finally,  if  we  are  willing  to  adjust 
the  model  parameter.-  rather  arbitrarily,  an  analytical  model 
can  be  found  that  will  have  nearly  the  same  equivalent 
source  as  that  for  the  finite  difference  earthquake. 

3 • 6 COMPARISON  WITH  DISLOCATION  MODELS 

For  the  relaxation  source  model  of  the  previous  sec- 
tion we  had  to  select  a set  of  model  parameters  that  seemed 
appropriate,  compute  the  radiation  field  and  then  see  how 
it  compared  with  our  source.  With  dislocation  models  we 
take  a different  approach.  From  the  radiation  field  we 
find  the  "equivalent  dislocation  model"  that  best  repre- 
sents the  finite  difference  earthquake.  The  parameters 
characterizing  this  equivalent  model  will  then  be  compared 

to  the  actual  parameters  in  the  finite  difference  calcula- 
tion. 


For  the  dislocation  model  geometry  we  select  the 
often  used  model  of  Haskell  (1964)  as  generalized  to  bi- 
lateral faulting  by  Savage  (1972).  if  we  consider  the  fault 


geometry  to  be  as  in  Figure  3.7,  a uniform  slip  is  assumed 
to  propagate  a distance  L in  the  +X_  direction  and  L in 

0 G IT 

the  -X„  direction  on  a rectangular  fault  of  length  L + L_ 

G o it 

and  width  w.  Then  the  amplitude  of  the  Fourier  transformed 
far-field  displacement  is 


|uC  (R,w)  I = M RC  ( R)  w |G(w)  I |FC(w,t  ,t  ) | , (3.11) 

***«■»  Q **  **  2 Q 7T 

where  c is  a or  8 corresponding  to  P or  S waves.  The  R 

has  as  components  the  spherical  coordinates  used  previously. 

Also, 

M = y w (L  + L ) D , (3.12) 

0 0*0 

where  D is  the  final  fault  offset  and 


Ra(R)  = e sin20  sin2<j>/  (4Trpa3R)  , 
*** 


(3.13) 


R6(R) 


[eQ  sin0  cos0  sin2<J>  + e^  sin0  cos2<j>]  / (4ttpB  3R)  , 


where  e , eQ,  e,  are  unit  vectors  in  the  R,0,<j>  directions 

— 2.  o 9 * 

and  p is  density.  The  G(w)  represents  the  Fourier  transformed 

rupture  time  history  and  F is  a factor  that  includes  the 

2 

modulation  of  the  radiated  waves  by  rupture  velocity  ef- 
fects. The  F is  given  by  Savage  (1972,  equation  10)  with 
2 

the  exception  that  his  cos0  should  be  replaced  by  sin0  cos<f> . 

An  equivalent  dislocation  model  for  the  finite  dif- 
ference earthquake  source  (FD)  should  have  the  same  moment 

(M  ) and  a time  history  of  faulting  (G(t))  that  gives  che 
o 

same  spectral  shape.  Let  us  look  first  at  the  momenc.  We 
previously  pointed  out  that  at  low  frequencies  the  FD  can 
be  represented  by  a double-couple.  In  fact,  for  w small 
(or  for  a pure  double-couple  source)  the  P and  SV  displace- 
ments are 
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-ik  R 


u (R,6,$,u)  = — § 
R k2  R 


B ^ 4 ^ (co ) P2  (cosQ)  sin24>  , 
2 2 2 


u fl(R,e,«D,w)  = — f 

o t2  R 

6 


-ikgR 


-A(1)  (u)+B(2)  (u) 
2 1 2 1 


(3.14) 


Pl  (cosQ)  sin<()Cos(|). 
2 


Note  that  for  a pure  double-couple  source  A ^ = - B^. 

2 12  1 

For  the  dislocation  model  Savage  (1972)  points  out 
that  approaches  unity  at  low  frequencies.  Then  equating 
(3.11)  and  (3.14),  for  small  frequencies  we  find  that 


M io)G(aj)  * 12irpa5  -3- 


B(4)  (w) 


00 


M iuG(ai)  = 12tt p 3 ' 
0 


-A(1)  (u>)  + B(2) 


JLL 


(a>) 


CO 


(3.15) 


The  moment  is  then  obtained  from  the  limit  as  oj  0.  Since 
G (co ) -*■  oj“l  and  A^  , B^  + go"2  (see  Figure  3.6)  in  the  low 
frequency  limit,  we  can  compute  the  moment  from  (3.15). 

From  the  P waves  we  find 


M = 1.62  * io 23  dyne-cm, 
0 

and  from  the  S waves 


(3.16a) 


M = 1.51  * 1023  dyne-cm.  (3.16b) 

0 

We  now  turn  our  attention  to  the  spectral  shape.  For 

convenience  we  normalize  the  amplitude  spectra  by  dividing 

by  the  moment  (M  ) and  the  double-couple  radiation  pattern 
c o 

(R  ) . The  normalized  far-f ield  spectra  for  the  dislocation 
model  are  then 
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(3.17) 


r 


i 


AC(0))  = 0)  J G (co ) I |f^  (oj,To  ,T^)  I , 


and  a similar  form  may  be  written  for  the  multipole  source 
representation.  Note  that  (3.17)  is  still  radiation  pattern 
dependent  via  which  plays  a role  similar  to  that  of  the 
higher  order  terms  (l  > 4)  in  the  multipolar  expansion. 

Two  forms  will  be  considered  for  the  rupture  time 
history.  First,  a form  suggested  by  Haskell  (1964) 


GH (t)  * 0,  t < 0 , 

GH(t)  = 1 - e‘t/T,  t > 0. 


(3.18) 


The  second  is  a slightly  more  complex  form  suggested  by 
Ohnaka  (1973) 


Gq (t)  =0,  t < 0 , 

Gq (t)  * 1 - (l-t/T)e"t/T,  t > 0. 


(3.19) 


If  the  rise  time  is  defined  as  the  time  to  reach  0.9,  we 
have  tR  = 2.3T  for  the  Haskell  form  and  tR  = 3.9T  for  the 
Ohnaka  form.  The  Fourier  transformed  amplitudes  are 


|Gh(oj)  I = w"1  (1  + u^T2)-1/2, 

A j 

|Gq  (<d)  | * CL)-  (1  + U)2T2)_1  . 


(3.20) 


If  we  select  the  elastic  velocities  (a, 8),  rupture 
velocity  (VR)  and  fault  length  (L  = = 0.5)  to  be  the 

same  as  in  the  finite  difference  calculation,  the  only  free 
parameter  for  matching  spectral  shapes  is  T.  For  each  of 
the  time  histories  of  (3.20)  a T was  selected  to  bring  the 
corner  frequencies  into  agreement,  using  the  half-amplitude 
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*:ule.  For  any  given  T the  agreement  is  better  in  some 
portions  of  the  radiation  pattern  than  in  others,  so  a 
value  that  "optimized"  the  agreement  over  the  pattern  was 
chosen.  It  was  previously  mentioned  that  the  ratio  of  P 
to  S wave  corner  frequencies  is  >1  for  our  finite  difference 
model  but  is  <1  for  dislocation  models  of  this  type.  To 
account  for  this  difference  a different  T was  selected  for 
the  P wave  spectra.  The  selected  "best"  values  for  the 
equivalent  'dislocation  are  as  follows:  For  the  S wave  and 

GH(t)  we  have  T = 0.19,  tR  = 0.44  and  for  GQ (t) , T = 0.11, 
tR  - 0.43.  For  the  P wave  and  G^(t)  we  have  T = 0.17,  t = 
0.39  and  for  GQ(t),  T = 0.10,  tR  = 0.39.  All  values  are^n 
seconds.  We  see  that  in  terms  of  rise  time  the  two  disloca- 
tion time  histories  give  the  same  result.  The  T values  are 
considered  to  be  + 0.01  for  a reasonable  fit  to  the  corner 
frequencies.  An  illustration  of  the  resulting  fit  is  shown 

in  Figure  3.12  for  S waves  at  one  location  on  the  radiation 
pattern. 

Having  specified  the  moment  and  rupture  time  history, 
we  have  an  equivalent  dislocation  model.  Let  us  now  compare 
the  parameters  of  this  model  with  the  actual  faulting  in  the 

difference  code.  First  consider  the  rise  time.  The 
shape  of  the  rupture  time  histories  shown  in  Figure  3.3  are 
compatible  with  the  dislocation  model  forms  of  (3.18)  or 
(3.19).  However,  the  rise  times  are  certainly  no  longer 
than  0.10  seconds  which  is  a factor  of  four  less  than  those 
for  the  equivalent  dislocation  model.  Can  this  discrepancy 
be  attributed  to  the  propagation  from  the  fault  surface  to 
the  elastic  radius?  The  displacement  time  histories  at  the 
elastic  radius  in  Figure  3.5  show  a rise  time  of  only  z 0.30 
seconds  even  though  the  interference  between  the  P and  S 
waves  seems  to  be  nearly  doubling  the  rise  time.  Another 
interesting  observation  is  that  the  rupture  time  (0.5/V  ) 
is  = 0.19  seconds,  a value  that  seems  consistent  with  the 
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Figure  3.12.  Comparison  of  the  far-field  SH  spectrum  at  t * 30°, 

4>  = 30° , for  the  finite  difference  earthquake  source 
to  spectra  from  a dislocation  model  with  two  rup- 
ture time  histories  (3.18,  3.19).  The  spectra  are 
normalized  by  the  moment  and  the  double-couple 
radiation  pattern  effect.  For  the  dislocation 
time  histories  T has  been  selected  to  bring  the 
comer  frequencies  into  agreement,  using  the 
half  amplitude  rule. 
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rise  time  in  Figure  3.5.  In  short,  we  are  forced  to  con- 
clude that  the  rise  time  for  our  equivalent  dislocation 
model  is  a factor  of  2-4  larger  than  the  rise  times  asso- 
ciated with  the  actual  rupture  process.  Perhaps  more  com- 
plex dislocation  model  geometries  and  rupture  time  histories 
would  reduce  the  discrepancy  but  removing  it  entirely  would 
appear  to  be  a formidable  task. 

What  about  the  seismic  moment?  We  obtained  two  esti- 
mates in  (3.16),  one  from  the  P waves  and  one  from  the  S 
waves.  The  dependence  of  moment  on  the  source  parameters 
is  given  by  (3.12).  If  we  take  the  fault  area  to  be 
0.6  x 1.0  km2  and  the  elastic  shear  modulus  to  be  p = 324 
kbar,  we  can  find  the  average  dislocation,  D , for  the 
equivalent  dislocation  model.  The  values  are  83.6  and  77.5 
cm  for  the  P and  S wave  moments  or  an  average  equivalent 
D 0 of  = 80  cm.  The  actual  displacement  at  five  positions 
on  the  fault  surface  was  shown  in  Figure  3.3.  Multiplying 
by  two  to  obtain  the  dislocation,  these  values  are  22-30 
cm.  In  fact,  the  average  dislocation  over  the  whole  fault 
surface  is  26  cm  as  is  shown  in  the  left  plot  in  Figure 
3.13.  The  actual  dislocation  is  then  smaller  than  the 
equivalent  dislocation  by  a factor  of  three. 

To  what  can  we  attribute  the  dislocation  discrepancy? 
We  can  hardly  suppose  that  the  effective  rigidity  in  the 
fault  zone  was  much  larger  than  the  e'.astic  value.  If  any- 
thing, the  plasticity  would  tend  to  have  the  opposite 
effect.  Another  possibility  is  that  the  equivalent  disloca- 
tion model  sees  a different  fault  area  than  that  on  which 
slip  occurs.  A nearly  maximum  effect  is  obtained  by  con- 
sidering the  effective  fault  surface  to  be  the  2.0  x 1.6  km2 
plane  one  zone  or  0.1  km  from  the  actual  fault  surface. 

With  this  fault  area  the  effective  dislocation  is  * 15  cm. 

The  actual  dislocation  on  this  plane  in  the  finite  difference 
grid  is  shown  in  the  right  of  Figure  3.13.  The  average  value 
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Figure  3.13.  The  static  dislocation  on  the  fault  surface  (left) 
and  one  zone  (0.1  km)  from  the  fault  surface.  The 
average  dislocation  (da)  over  the  plane  is  indi- 
cated by  a dotted  line  on  each  plot.  The  data 
points  are  the  values  at  the  nodes  in  the  grid. 

The  nodes  are  indicated  by  their  distance  from 
the  X = 0 or  Z = 0 line. 
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is  10.6  cm  and  the  discrepancy  has  been  reduced  to  a more 
satisfactory  factor  of  1.4.  Some  further  improvement 
could  probably  be  made  by  taking  an  optimum  portion  of  this 
plane.  However,  we  must  address  the  question:  When  we 

construct  an  equivalent  dislocation  model  for  an  earthquake, 
what  do  we  mean,  in  terms  of  the  physics,  by  the  fault  sur- 
face and  average  dislocation? 

The  product  of  aiea  times  average  static  displacement 
for  the  two  planes  in  Figure  3.13  is  greater  for  the  plane 
away  from  the  fault  surface.  This  kind  of  radiation  of  dis- 
placements from  the  fault  surface  is  unlike  that  in  a dis- 
location model  and  seems  to  be  due  to  the  plastic  adjust- 
ment of  the  stress  concentrations  that  arise  near  the  rup- 
ture front.  More  calculations  must  be  done  to  be  sure  we 
understand  the  effect  of  finite  difference  zone  size,  details 
of  the  material  behavior  at  high  stress,  etc.,  but  our  model 
of  the  decay  of  displacements  from  the  fault  surface  seems 
more  likely  to  correspond  to  physical  reality  than  that  in 
the  dislocation  model. 

3.7  CONCLUSIONS 

Our  purpose  in  undertaking  this  work  was  twofold. 
First,  we  wanted  to  carry  out  a three-dimensional  simulation 
of  earthquake  faulting  and  to  link  the  finite  difference 
code  with  elastic  wave  theory;  in  other  words,  to  explore 
the  feasibility  of  the  technique.  Second,  we  wanted  to  see 
what  could  be  learned  about  earthquake  source  theory  from 
the  results  of  our  first  calculation.  Our  concluding  re- 
marks are  directed,  in  turn,  to  these  two  purposes. 

The  event  selected  for  study  was  a small  bilateral 
earthquake  with  a symmetrically  propagating,  constant 
velocity  rupture.  A simple  mechanism  for  relieving  the 
stress  concentrations  at  the  crack  tip  was  introduced  by 
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f considering  the  material  in  the  fault  zone  to  behave  in  an 

elastic-plastic  manner.  What  seems  to  be  a realistic  simu- 
lation of  stick-slip  fault  rupture  was  the  result.  This 
numerical  simulation  is  perhaps  one  of  the  simplest  that 
can  be  done.  The  point  is  that  now  much  more  complex  geo- 
metries and  material  behaviors  can  be  introduced  and  their 
effects  determined. 

propagating  the  earthquake  generated  waves  very 
far  from  the  fault,  much  efficiency  is  gained  by  represent- 
ing the  event  by  an  equivalent  elastic  source.  A generally 
applicable  technique  for  obtaining  such  a source  from 
finite  difference  output  was  developed  during  the  course  of 
this  research.  The  equivalent  elastic  source  representation 
of  finite  difference  source  calculations  also  has  great 
value  in  allowing  convenient  comparison  of  the  sourne 
characteristics  with  those  obtained  from  other  theories. 

We  have  presented  the  results  of  a single  calculation; 
there  are  no  others.  With  due  caution  against  pushing  the 
results  too  far,  we  find  a number  of  intriguing  results 
when  comparing  the  far-field  displacement  spectra  from  our 
finite  difference  source  to  that  from  the  eiastodynamic 
model  of  Archambeau/Minster  (Archambeau,  1968;  Minster,  1973) 
and  to  that  from  dislocation  models  (Haskell,  1964;  Savage, 
1972) . In  almost  all  respects  the  radiation  from  our 
source  more  closely  resembles  that  from  the  Archambeau/ 
Minster  volume  source,  relaxation  model.  The  long  period 
level  or  moment  is  essentially  identical.  At  frequencies 
up  to  the  corner  frequency  the  radiation  patterns  are  also 
nearly  the  same.  While  the  corner  frequencies  themselves 
and  the  shape  of  the  spectra  in  this  region  do  differ  for 
the  analytical  a)  id  finite  difference  sources,  it  is  interest- 
ing to  note  tha:  the  ratio  of  P to  S corner  frequencies  is 
almost  the  same  over  the  whole  radiation  pattern  and  is 
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greater  than  one  everywhere.  In  short,  the  elastodynamic 
relaxation  source,  using  the  same  model  parameters  as  in 
our  finite  difference  source,  is  quite  similar  as  a radiator 
of  far-field  seismic  waves. 

Using  the  dislocation  theory  of  Haskell  (1964)  ex- 
tended to  bilateral  faulting  by  Savage  (1972) , we  constructed 
an  equivalent  dislocation  model  from  the  far-field  radiation 
computed  for  our  fin: te  difference  source.  Optimizing  the 
agreement  between  corner  frequencies,  we  found  that  the  rise 
time  for  the  dislocation  model  time  history  was  a factor  of 
2-4  greater  than  the  actual  rise  time  for  the  displacement 
time  histories  on  the  fault  surface.  Examination  of  the  time 
histories  from  which  the  equivalent  elastic  source  was  ob- 
tained indicated  that  this  was  not  due  to  finite  grid  attenua- 
tion or  dispersion  effects.  Using  the  actual  fault  surface 
and  elastic  shear  modulus,  the  average  fault  offset  for  the 
equivalent  dislocation  model  was  more  than  a factor  of  3 
greater  than  the  fault  offset  in  ttie  finite  difference  cal- 
culation. 

In  terms  of  a kinematic  description  of  the  displace- 
ments at  the  fault  surface,  our  finite  difference  earthquake 
looks  quite  similar  to  a dislocation  model.  Nevertheless, 
we  found  rather  poor  agreement  between  the  actual  fault  para- 
meters and  those  of  the  equivalent  dislocation  model  when 
we  consider  the  absence  of  errors  due  to  propagation  and 
sampling  effects. 

Why  does  the  finite  difference  source  look  more  like 
a volume  source  model?  It  is  tempting  to  think  that  it  is 
because  the  entire  "fault  zone"  is  yielding  in  the  same  way 
as  the  spherical  volume  in  the  Archambeau/Minster  source. 

But  this  seems  to  be  not  the  case.  Both  P and  S waves  are 
radiated  directly  from  the  fault  surface  (Figure  3.5).  The 
plastic  work  (stress  time  plastic  strain)  in  the  fault  zone 
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is  large  in  the  vicinity  of  the  crack  tip  during  rupture. 

But  in  most  of  the  fault  zone  the  amount  of  plastic  flow  is 
guite  small.  We  saw  in  Figure  3.13  that  the  plastic  re- 
adjustment of  the  stress  concentrations  near  the  rupture 
front  cause  the  displacements  one  zone  from  the  fault  sur- 
face to  be  quite  different  than  they  would  be  in  an  elasto- 
dynamic  dislocation  model.  This  seems  to  be  the  cause  of  the 
disagreement  between  equivalent  dislocation  and  actual  model 
parameters.  However,  it  is  not  entirely  clear  why  the 
volume  source  model  agrees  so  well.  Nor  do  we  know  whether 
this  agreement  is  particular  to  this  finite  difference 
earthquake  or  would  be  obcained  for  a whole  class  of  finite 
difference  earthquake  simulations.  Answers  to  these  ques- 
tions will  emerge  from  the  study  of  future  calculations. 


IV.  SHORT  PERIOD  THEORETICAL  SEISMOGRAMS  FOR  A 
PROPAGATING  SOURCE  MODEL  FOR  EARTHQUAKES 

4.1  INTRODUCTION 

Much  of  the  available  evidence  for  the  character  of 
the  earthquake  source  mechanism  is  contained  in  recordings 
of  ground  motion  at  great  distances  from  the  source.  The 
same  can  be  said  for  underground  nuclear  explosions.  The 
far-field  seismic  recordings  are  often  quite  complex  with 
the  causes  lying  in  two  categories:  (1)  Complexity  of  the 

source/  and  (2)  Modulation  of  the  source  generated  ground 
motion  by  the  earth  structure  intervening  between  source 
and  receiver.  For  nearly  all  seismological  studies  we  want 
to  be  able  to  remove  one  of  these  causes  cf  complexity  in 
order  to  understand  the  other. 

In  a previous  paper,  Bache  and  Harkrider  (1976,  sub- 
sequently referred  to  as  paper  I)  set  forth  a method  by 
which  high  resolution  methods  for  simulating  the  complex 
seismic  source  can  be  linked  with  high  resolution  elastic 
wave  propagation  techniques.  The  particular  application  ad- 
dressed in  that  paper  was  to  the  body  waves  propagating  to 
the  far-field.  For  surface  waves  a similar  formulation  was 
given  by  Harkrider  and  Archambeau  (1976). 

In  this  paper  the  objective  is  to  use  the  theoretical 
results  given  in  paper  I to  compute  synthetic  seismograms 
for  an  earthquake  source  model  that  is  sufficiently  complex 
to  illustrate  the  potential  utility  of  this  theory  for 
seismological  research.  The  earthquake  source  model  chosen 
is  that  of  Archambeau  (1968)  and  Minster  (1973).  Most 
prominent  among  the  "higher  order"  effects  included  in  the 
version  of  the  model  used  here  are  those  due  to  fault  propa- 
gation at  finite  rupture  velocity  and  these  can  have  important 
effects  on  the  teleseismic  short  period  ground  motion  record- 
ings. 
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We  begin  by  describing  the  rapid  and  inexpensive 
techniques  being  used  for  computing  synthetic  seismograms. 
These  techniques  include  the  interfacing  of  the  paper  I 
formulation  for  the  steeply  emergent  body  waves  at  the  base 
of  the  crustal  model  at  the  source  with  generalized  ray 
theory  (e.g.,  Wiggins  and  Helmberger,  1974)  for  the  upper 
mantle  travel  path. 

The  Archambeau/Minster  source  model  is  briefly  dis- 
cussed/ particularly  the  scaling  with  the  source  parameters. 
Our  intention  in  this  section  is  p,  imarily  to  demonstrate 
the  kinds  of  effects  that  can  now  be  conveniently  included 
in  theoretical  seismograms  rather  than  to  present  a detailed 
study  of  these  effects.  Therefore/  only  a single  source  at 
one  orientation  (strike-slip)  is  used  for  the  theoretical 
seismogram  calculations.  For  this  source  the  dependence  of 
teleseismic  body  wave  recordings  on  focal  depth  and  fault 
azimuth  is  discussed. 


4.2  COMPUTATION  OF  THEORETICAL  SEISMOGRAMS 


The  Fourier  transformed  equations  of  motion  in  a 
homogeneous,  isotropic  elastic  medium  may  be  written 


where  u is  particle  displacement  and  k and  k„  are 
, a B 

the  compressional  and  shear  wave  numbers.  The  Cartesian 

_(4)  _ 

potentials  x and  x can  be  written  in  terms  of  the 
following  expansion  in  spherical  harmonics; 


00  I 

X^(R,w)  = h^2^  (kjR)  f A^  (to)  cos  m^i 

£=0  m=0 

+ Bim  (u>)  sin  m<*]  pIJ(c°s0),  j = 1,2, 3,4  , (4.2) 
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where  the  h^  are  spherical  Hankel  functions  of  the 


second  kind  and  the  P 


.m 


^ are  associated  Legendre  functions. 
The  vector  R has  as  components  the  spherical  coordinates 
R,0,4> 


As  was  pointed  out  in  paper  I , nearly  any  proposed 
seismic  source  model  can  be  specified  in  terms  of  the  ex- 
pansion {multipole)  coefficients  of  (4.2).  Therefore,  (4.2) 
together  with  (4.1)  provide  an  elastic  point  source  repre- 
sentation of  the  (outgoing)  displacement  field.  For  a source 
prescribed  in  this  form  and  buried  at  a depth  h in  a plane- 
layered, serai-infinite  elastic  medium,  paper  I includes  the 
derivation  of  the  expressions  for  the  steeply  emergent  body 
waves  emanating  into  the  underlying  halfspace  at  a specified 
angle  0n  or,  equivalently,  horizontal  phase  velocity  c. 

The  geometry  is  shown  in  Figure  4.1.  Polarized  into  P,  SV 
and  SH  components,  these  displacements  are  fEo.  (33), 
paper  I) 
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where  all  ..otation  is  defined  in  paper  I. 


To  compute  theoretical  seismograms  we  take  the  pi 
layered  model  of  Figure  4.1  to  represent  the  crust  in  the 
source  vicinity.  Then  (4.3)  gives  the  waves  propagating 
into  the  upper  mantle.  Other  methods  are  then  used  to 
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Figure  4.1.  The  geometry  and  coordinate  system  for  a 
at  depth  h ir  a multilayered  half space 
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source 
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compute  the  remainder  of  the  travel  path.  While  there  are 
alternatives,  we  have  found  it  convenient  to  use  generalized 
ray  theory  (Wiggins  and  Helmberger,  1974)  for  the  upper 
mantle  and  Haskell-Thomson  matrices  (Haskell,  1960,  1962) 
for  the  crust  in  the  receiver  vicinity.  Thus,  the  travel 
path  is  divided  into  three  pieces: 

The  crust  at  the  source  to  a depth  D. 

The  crust  at  the  receiver  to  the  same  depth  D. 

A laterally  homogeneous  upper  mantle  extending 
from  D to  the  deepest  depth  of  interest. 

The  travel  path  segments  are  linked  by  requiring  that 
the  velocities  at  D be  the  same  in  all  three  structures. 
This  scheme  gives  a great  deal  of  computational  flexibility 
in  that  portions  of  the  travel  path  can  be  varied  without 

repeating  the  entire  calculation.  The  actual  implementation 
is  as  follows: 


1. 

2. 

3. 


A generalized  ray  theory  calculation  is  done  to 
determine  the  impulse  response  at  epicentral  dis- 
tances of  interest.  The  source  and  receiver  are 
at  the  surface  of  an  earth  model  in  which  the  top 
D kilometers  have  been  removed.  A ray  parameter 
or  horizontal  phase  velocity  (c)  is  associated 
with  each  major  arrival  at  the  receiver  and  this 
is  an  input  parameter  "nr  the  crustal  response 
calculations,  in  practic  the  crustal  response 
variation  with  c is  sufficiently  small  that  for 
most  distances,  certainly  all  those  beyond  the 
triplications,  a single  c is  sufficient. 

The  displacements  exiting  the  source  region  at 
phase  velocity  c are  computed  using  (4.3). 

A frequency  dependent  transfer  function  is  com- 
puted for  the  plane  layered  model  of  the  crust 
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beneath  the  receiver.  The  appropriate  transfer 
functions  were  given  by  Haskell  (1960,  1962). 
For  example,  for  P waves  the  transfer  function 
is 


Ty  (u>)  = 2c 


(jr  - jr  ) 

\ 4 1 3 1 / 

a D_ 
n J 


for  the  vertical  component  of  motion  at  the  surface 
and 


Tr(“) 


2c 


a 

n 


) 


$ 


for  the  radial  component  (all  notations  are  as  in 
paper  I) . Here 


Anelastic  attenuation  and  the  attendant  dispersion 
are  included  via  the  operator 


exp 


Jin 


(X000)]) 


f 


using  the  formulation  suggested  by  Strick  (1970) . 
Here  f is  frequency  and  the  controlling  parameter 
is  T/Q,  the  ratio  of  travel  time  to  the  average 
path  material  quality  factor.  The  T/Q  is  conven- 
iently input  as  a constant  prescribed  for  each 
seismogram  calculation.  However,  it  can  easily 
be  computed  from  a postulated  Q versus  depth 
model  for  each  ray  path.  Also,  Hasegawa  (1971) 
has  pointed  out  that  only  minor  modifications  are 
required  to  include  the  effect  of  Q in  Haskell- 
Thomson  matrix  calculations. 
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The  sensing  instrument  is  characterized  by  a fre- 
quency dependent  transfer  function. 

The  source,  source  crust,  receiver  crust,  Q and 
instrument  effects  are  computed  in  the  frequency  domain. 
Numerical  Fourier  transformation  brings  their  combined  ef- 
fect to  the  time  domain  where  it  is  convolved  with  the 

upper  mantle  response  to  complete  the  theoretical  seismo- 
gram calculation. 

In  later  sections  we  will  present  theoretical  seismo- 
grams computed  with  this  method.  However,  we  first  describe 
the  source  model  and  its  radiation  field. 

4-3  THE  ARCHAMBEAU /MINSTER  EARTHQUAKE  SOURCE  Hnnyr. 

Archambeau  (1968)  and  Minster  (1973)  have  developed  a 
relaxation  source  theory  for  earthquakes  in  which  the  finite 
rupture  velocity  causes  significant  excitation  of  higher 
order  (I  > 2)  terms  in  the  multipolar  expansion  (4.2)  of  the 
radiation  field,  particularly  at  frequencies  beyond  the 
corner  frequency.  These  higher  order  terms  have  an  inter- 
esting effect  on  the  body  waves  propagated  to  teleseismic 
distances. 

The  main  features  of  the  elastic  radiation  field  from 
this  earthquake  source  model  are  discussed  by  Minster  and 
Archambeau  (1976) . For  our  calculations  we  shall  use  the 
model  they  call  the  "tangentially  expanding  sphere”.  In 
this  model  a spherical  volume  of  vanishing  rigidity  expands 
rom  zero  to  its  final  diameter  with  a point  on  the  diameter 
ixed  in  space.  This  expanding  weak  zone  is  assumed  to  be 
m a homogeneous  elastic  medium  that  is  stressed  in  pure 
shear.  The  parameters  of  the  earthquake  model  are  then: 

L “ The  fault  length:  equal  to  the  final 
diameter  of  the  expanding  sphere. 
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VR  - The  rupture  velocity;  controls  the  rate 
of  increase  of  the  sphere. 

Qt/B  — The  P and  S elastic  wave  velocities. 

S {0)  - The  prestrain;  equal  to  a^°Vu>/  the  ratio 
of  prestress  to  material  rigidity.  The 
radiation  field  scales  directly  with  S^. 

Rg  — The  radius  of  the  zone  in  which  the  tectonic 
prestress  is  released. 

Archambeau's  assumption  of  a finite  Rg,  which  is  an 
approximation  met nt  to  represent  a nonuniform  prestress 
field,  has  led  to  considerable  controversy  (e.g.#  Molnar, 
et  al . / 1973)  since  the  resulting  source  spectrum  is  peaked. 
While  the  assumed  value  of  Rg  has  great  importance  for 
long  period  studies  or  very  large  faults,  its  value  has 
almost  no  effect  on  short  period  body  waves  for  moderate 
fault  dimensions.  Here  we  shall  simply  take  Rc  to  bo 
infinite. 

While  the  scaling  of  the  radiated  field  with  these 
parameters  is  discussed  in  detail  in  Minster  and  Archambeau 
(1976),  we  should  point  out  several  features  that  will  be 
of  importance  here.  For  a particular  fault  orientation, 
the  multipole  coefficients  are  completely  specified  by 
Poisson's  ratio  (v) , prestrain  and  the  ratios  L/V_  and 

X\ 

Vr/B  (assuming  Rg  ■+■  »)  . The  amplitude  of  the  coefficients 
is  fixed  by  L/B  (or  L/a  since  v is  fixed)  while  the 
corner  frequency  depends  on  L/V R.  The  character  of  the 
radiated  field  also  depends  on  the  orientation  of  the  pre- 
strain and  rupture  direction  with  respect  to  the  free  sur- 
face. These  can  conveniently  be  specified  by  strike,  dip 
and  plunge  angles. 

The  amplitude  dependence  on  L/8  is  identical  to 
the  behavior  of  elementary  dislocation  theories.  For  example. 
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given  in  paper  I are  the  multipole  coefficients  for  a hori- 
zontal double-ccuple  (normal  strike-slip  fault)  dislocation 
source.  A typical  term  is 


a'1’ 


M 


i D (w)  kjjAf 
24i 


(4.4) 


where  D(w)  is  the  dislocation  per  unit  area  and  A^  is 
the  fault  area.  If  the  fault  dimension  is  » L and  the 
strain  drop  is  approximately  D/L,  (4.4)  shows  that 
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(4.5) 


which  is  the  same  amplitude  scaling  that  occurs  in  the 
Archambeau/Minster  model. 

We  shall  compute  theoretical  seismograms  for  a parti- 
cular. earthquake  specified  by  the  following  source  parameters: 

VR  = 4.14  km/sec, 

6 = 4.6  km/ sec,  (4.6) 

RS  " "• 

The  fault  orientation  is  assumed  to  be  right— lateral  strike- 
slip  along  the  geographic  X^-axis,  assumed  to  point  north. 

The  geographic  c" ordinate  system  and  fault  geometry  are  shown 
in  Figure  4.2.  Also  shown  are  the  azimuth  and  takeoff  angle 
(<pf  t)  foi  specifying  a ray  (A).  In  terms  of  the  coordinates 
of  (4.2)  and  Figure  4.1,  t may  be  identified  with  6g 
measured  from  the  downward  vertical  in  the  source  layer. 

Then  t = 0g  = tt  - 0 where  9 is  the  spherical  coordinate 
in  (4.2). 

In  deriving  the  steeply  emergent  body  waves  given  by 
(4.3),  only  the  far-field  terms  (decaying  as  R-1)  were  re- 
tained. In  the  far-field  the  dependence  of  the  expansion 
(4.2)  on  R is  separable  since  the  Hankel  functions  reduce 


L = 10  km. 


a = 7.9  km/sec. 


(0) 


10 


-3 
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(Up) 


Fault  Plane 


Figure  4.2.  The  geometry  and  coordinate  system  for  a right- 
lateral  strike -slip  fault.  The  outgoing  dis- 
placement field  along  the  ray  A is  specified 
by  the  azimuth  <J>  and  takeoff  angle  x. 
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1 ft-  R » 3 = ci  or  0. 


(4.7) 


Further,  the  P and  S waves  separate  with  the  P waves  given 
by  tae  first  term  on  the  right  side  of  (4.1)  and  the  S waves 
given  by  the  second  term. 


In  view  of  the  similarity  properties  of  the  model,  the 
multipole  coefficients  for  the  earthquake  specified  by  (4.6) 
are  identical  to  those  for  an  event  with  fault  parameters: 


L = 7.6  km. 


a = 6.0  km/ sec. 


S 


(0) 


VR  = 3.14  km/sec, 
8 = 3.49  km/sec. 


(4.8) 


While  the ^multipole  coefficients  are  identical,  the  stress 
drop  (p  S ) differs  depending  on  the  shear  moduli.  From 
(4.1),  it  is  clear  that  tha  far-field  P and  S wave  displace- 
ments are  proportional  to  a2  and  32,  respectively. 

Before  presenting  theoretical  seismograms  for  the 
strike-slip  earthquake  source  (4.6)  or  (4.8),  it  is  useful 
to  develop  some  insight  into  the  behavior  of  the  source  as 
a radiator  of  elastic  waves  in  a whole  space.  Far-field 
displacement  radiation  patterns  are  presented  in  Figures 
4.3  and  4.4.  These  are  computed  from  (4.1)  and  (4.2),  re- 
taining only  far-field  terms,  using  the  source  described 
by  (4.6).  The  patterns  are  equally  valid  for  the  source 
(4.8)  with  the  amplitudes  being  scaled  by  the  square  of 
the  velocity  ratios  or  (6. 0/7. 9) 2 . 

In  Figure  4.3  the  far-field  displacement  radiation 
patterns  are  shown  for  a number  of  frequencies  and  several 
levels  of  truncation  for  the  infinite  series  (4.2).  These 
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Figure  4.3.  Far-field  displacement  radiation  patterns  at  three 
frequencies  (plotted  vertically)  and  three  levels 
truncation  of  the  multipolar  coefficients.  The 
patterns  are  for  fixed  azimuth  <j>  = 45°  and  are 
oriented  with  respect  to  the  system  of  Figure  4.2 
as  indicated  in  the  double-couple  patterns  of  the 
first  row.  The  SH  component  is  suppressed  on  the 
N - 4 and  first  set  of  N * 6 plots.  (The  meaning 
of  the  plot  symbols  is  indicated  on  the  last  two 
plots  of  the  center  row.)  The  displacement  ampli- 
tude at  the  outer  circle  is  indicated  in  units  of 
a,  a = 1.98  cm-sec. 


radiation  patterns  are  for  the  three  components  (u^  or  P, 

Uq  or  SV  and  u0  or  SH)  of  the  amplitude  j u(R,  0 , <J>  ,w)  | of 
the  displacement  spectrum  at  fixed  R,4>,u)  and  varying  0. 

The  azimuth  is  <p  = 4 5°  and  we  have  chosen  the  frequencies 
f = 0.1,  0.5,  1.0  Hz.  The  spectral  amplitudes  are  for 
R = 100  km,  though  in  view  of  (4.7)  they  can  easily  be 
scaled  to  any  R of  interest.  Truncating  the  series  (4.2) 
at  £ = 2,  the  source  is  a quadrupole  or  double-couple.  We 
then  show  the  effect  of  adding  the  higher  order  terms  as 
we  truncate  at  £ - 4 and  l = 6.  The  level  of  truncation  or 
maximum  £ is  denoted  N.  For  the  double-couple  (N  = 2) 
the  <p  = 45°  orientation  is  a node  on  the  SH  wave  pattern. 
However,  a strong  SH  component  appears  upon  adding  higher 
order  terms.  For  clarity  this  component  is  suppressed  on 
the  N * 4 patterns  and  the  N = 6 patterns  are  shown  with 
and  without  the  SH  waves. 

From  Figure  4.3  we  see  that  at  low  frequencies  the 
source  is  rather  well  approximated  by  a double-couple.  At 
higher  frequencies  the  higher  order  terms  have  a pronounced 
effect,  skewing  the  radiation  pattern  toward  the  direction 
of  fault  propagation.  The  large  SH  component  is,  of  course, 
oniy  present  due  to  the  higher  order  multipoles.  The  0.5 
and  1.0  Hz  patterns  are  typical  of  those  controlling  tele- 
seismic  bodywaves.  Taking  higher  order  propagation  effects 
into  account,  these  patterns  are  certainly  quite  different 
from  those  fcr  a pure  double-couple. 

Further  insight  into  the  properties  of  the  source 
as  an  elastic  wave  radiator  is  provided  by  the  spectra  of 
Figure  4.4.  For  the  three  displacement  components  (P,  SV, 
SH)  spectra  are  shown  at  a forward  (<p  = 45°)  and  backward 
(<P  = 135°)  azimuth  for  the  takeoff  angle  x = 30°.  For 
these  spectra  the  series  truncation  is  at  N = 6.  The 
measured  (from  the  high  and  low  frequency  asymptotes)  P 


and  S corner  frequencies  for  the  <p  = 45°  spectra  are 

f0(P)  * °*3  Hz  and  f0<S)  * 0-22  Hz.  We  mention  in  passing 
that  these  are  in  rough  agreement  with  empirically  deter- 
mined values  for  actual  earthquakes  with  similar  source 
dimensions  {e.g.,  Wyss  and  Shamey,  1975). 


Below  the  corner  frequency,  the  source  is  pre- 
dominantly double-couple.  At  greater  frequencies  the 
higher  order  terras  assume  increasing  importance  and  come 
to  dominate,  as  is  apparent  in  the  radiation  patterns.  In 
fact,  truncating  the  series  after  six  terms,  convergence 
has  not  yet  been  attained  for  frequencies  much  beyond  1 Hz. 
Minster  (private  communication)  has  found  that  the  addition 
of  more  terras  acts  to  smooth  the  spectrum,  removing  the 
spectral  holes.  However,  the  theoretical  seismograms  to  be 
presented  are  insensitive  to  details  of  the  source  spectrum 

beyond  1 Hz  and  we  need  not  be  concerned  with  the  complicated 
structure  in  this  spectral  region. 


4 • 4 THEORETICAL  SEISMOGRAMS 

We  now  specify  an  earth  model  and  compute  theoretical 
teleseismic  short  period  seismograms  for  thi  earthquake 
source  described  in  the  previous  section.  The  source  crustal 
model  is  representative  of  the  Basin  and  Range  Province  in 
the  western  United  States  (Hill  and  Pakiser,  1967)  and  is 
tabulated  m Table  4.1.  The  crust  at  the  receiver  (tabu- 
lated in  Table  4.2)  is  an  average  continental  crust  which 
has  almost  no  effect  on  the  seismogram  shape. 

The  upper  mantle  model  is  HWNE  (Helmberger  and 
Wiggins,  1971)  with  the  top  50  km  removed  for  compatibility 
with  our  computational  scheme.  The  generalised  ray  theory 
calculation  was  done  for  a range  of  4000  km.  For  the  total 
epicentral  distance  it  is  neiessary  to  add  corrections  for 
the  travel  in  the  crust  at  the  two  ends.  At  this  distance 
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Table  4.1.  Basin  and  Range  Crustal  Structure 


Depth 

a 

8 

P 

y 

(km) 

(km/sec) 

(km/sec) 

(gm/cm3 ) 

(kbar ) 

V 

2.5 

4.2 

2.6 

2.40 

162 

0.189 

18.0 

6.  0 

3.49 

2.82 

343 

0.244 

31.0 

6.  / 

3.85 

3.00 

489 

0.254 

46.0 

7.9 

4.6 

3.30 

698 

0.244 

50.0 

8.0 

4.6 

3.30 

698 

0.253 

Table  4.2.  Average  Continental  Crust 


Depth 

a 

8 

o 

(km) 

(km/sec) 

(km/sec) 

(gm/cm3) 

22.0 

6.03 

3.53 

2.8 

37.0 

6.70 

3.80 

3.0 

50.0 

8.00 

4.60 

3.3 

the  dominant  energy  arrives  with  horizontal  phase  velocity 
c = 12.99  km/sec,  the  value  used  for  the  crustal  structure 
calculations . With  this  c the  distance  correction  for  the 
receiver  crust  is  s 30.7  km  while  the  source  crust  correction 
is  generally  smaller,  depending  on  source  depth.  At  4000  km 
the  upper  mantle  response  is  little  different  from  a constant 
spreading  factor  with  a value  of  about  0.95  x 10~\  Finally, 
the  anelastic  attenuation  operator  was  applied  with  T/Q  = 

1.05  and  the  ground  motion  was  filtered  by  the  response  of 
a standard  LRSM  short  period  Benioff  seismometer. 
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Setting  up  the  theoretical  calculations  this  way, 
only  the  source  itself  and  its  modulation  by  the  local  geo- 
logy varies  from  one  theoretical  seismogram  to  another. 

That  is,  all  features  of  the  computations  are  held  fixed 
except  the  ground  motion  propagating  from  the  base  of  the 
crust  at  the  source  which  is  computed  from  (4.3).  For 
the  earthquake  source  specified  by  the  source  parameters 
(4.6)  or  (4.8)  and  the  orientation  of  Figure  4.2,  we  are 
interested  in  the  dependence  of  teleseismic  ground  motion 
on  azimuth,  focal  depth  and  the  proximity  of  major  velocity 
discontinuities  in  the  crust.  While  only  the  strike-slip 
orientation  is  being  considered,  we  should  point  out  that 
one  of  the  advantages  of  the  multipolar  source  representa- 
tion is  the  availability  of  computational  algorithms  (e.g.. 
Minster,  1973)  for  conveniently  rotating  the  source  to  any 
desired  orientation. 

We  now  consider  the  earthquake  source  at  focal  depths 
of  10,  15,  25,  35,  45  km  in  the  structure  of  Table  4.1.  For 
each  depth  it  is  the  multipole  coefficients  that  remain 
constant  and  the  source  parameters  then  change  according  to 
the  similarity  properties  of  the  model.  The  appropriate 
parameters  are  (4.6)  and  (4.8)  and  one  intermediate  to  these 
and  are  summarized  in  Table  4.3.  The  stress  drop  for 
S<°>  = 10  is  given  as  well  as  the  strain  drop  correspond- 
ing to  events  with  constant  stress  drop,  cr  ^ = 100  bars. 
Since  ;ach  seismogram  scales  directly  with  S ^ or  it 

is  easy  to  correct  to  other  convenient  values.  The  radia- 
tion patterns  and  spectra  of  Figures  4.3  and  4.4  apply 
equally  well  to  any  of  these  events,  being  scaled  by  the 
square  of  the  velocity  ratios.  Note  that  for  fixed  strain 
drop,  this  is  little  different  from  scaling  by  stress  drop 
ratios  since  = pB2  S^. 

The  synthetic  seismograms  are  shown  in  Figures  4.5 
and  4.6.  These  and  all  subsequent  seismograms  are  for  a 


Table  4.3. 


Source  Parameters  for 


Events  at  Five  Focal  Depths 


Source 
Depth 
h (km) 

Fault 
Length 
L (km) 

Rupture 
Velocity 
VR (km/sec) 

Stress  Drop 
a(0)  (bars) 

Strain 
a<°>  = 

Drop  for 
100  bars 

10,  15 

7.6 

3.15 

343 

0.29 

x 10~3 

25 

8.4 

3.47 

489 

0.20 

x 10“3 

35,  45 

10.0 

4.14 

698 

0.14 

x 10~3 

constant  stress  drop  of  100  bars.  In  Figure  4.5  are  shown 
the  records  for  events  at  the  five  focal  depths  and  at  two 
azimuths,  0 = 45“  and  135  . The  focal  depth  and  fault  length 
are  indicated  on  each  record  together  with  the  the  phase 
from  which  the  amplitude  measurement  was  made  and  the  ap- 
parent period  of  this  phase.  We  have  attempted,  as  much  as 
possible,  to  analyze  these  records  following  standard  field 
practice.  The  was  computed  from  m^  = log  A/T  +3.3 
where  A is  measured  peak-to-peak.  The  major  influerce  of 
the  apparent  period,  T,  is  in  correcting  to  "true"  ground 
motion  since  the  instrument  correction  is  large  and  rapidly 
varying  in  this  period  range. 

In  Figure  4.S  the  azimuthal  variation  is  shown  for 
the  25  km  focal  depth  event.  Only  azimuths  from  0°  to  180° 
are  shown  since  the  pattern  is  symmetric  about  the  north- 
south  (XG)  axis  except  for  a change  in  sign.  Of  particular 
interest  are  the  seismograms  near  <p  = 90°  where  the  sense 
of  first  motion  reverses.  This  azimuth  is  a node  on  the 
pattern  for  a simple  double-couple  representation  of  a 
strike-slip  fault. 


b.  Azimuth,  <p  = 135° 


ofe°rIn^raLSeiJra?9raraS.at  an  ePicentral  distance 
?fr.\4051G.kin  and  tw0  azimuths  with  respect  to  the 
strike-slip  earthquake  sources  of  Table  4.3  All 
records  have  been  scaled  to  stress  drop  o(0)  = iqc 
bars.  The  focal  depth,  fault  length  and  mb  are 

pefk-?oeLak  ^C?^e5°rd-  At  left  is  the  *aximum 
£!la\t0  ?ea^  amplitude  in  millimicrons  at  1 Hz 

mentis  ^ndicates  the  phase  at  which  the  mb  me,’sur 

ade  and  T ls  the  apparent  period  of  this 
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Theoretical  seismograms  at  15  azimuths  (0  < $ < 180) 
for  the  L = 8.4  km  source  at  a 25  km  focal  depth. 

The  m^  phase  from  which  it  was  computed,  apparent 
period  of  this  phase  and  maximum  peak-to-peak  ampli- 
tude (1  Hz)  are  shown  as  in  Figure  4.5. 
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MALYSIS  OF  THEORETICAL  SEISMnftPUMc 


Our  intention  here  is  primarily  to  demonstrate  the 
kinds  of  effects  that  can  be  included  in  the  theoretical 
seismogram  calculations.  A detailed  study  of  these  effects 
and  their  implications  for  source  mechanism  studies  is  now 
underway.  However,  several  interesting  aspects  of  the  theo- 
retical seismograms  of  Figures  4.5  and  4.6  will  be  discussed. 
In  particular,  we  address  the  questions:  (1)  What  is  the 

influence  of  the  propagation  effects  introduced  by  the  higher 
order  multipoles?;  and  (2)  How  are  the  seismograms  influenced 
y details  of  the  crustal  structure  in  the  source  region? 

In  Figure  4.7  we  compare  theoretical  seismograms  which 
differ  only  in  the  level  at  which  the  multipolar  expansion  is 
truncated.  The  N = 6 seismogram  also  appears  in  Figures  4.5 
and  4.6  where  it  is  marked  with  an  asterisk.  At  this  azimuth 
the  most  important  effect  of  the  higher  order  terms  seems  to 
..a  on  the  frequency  content  of  the  seismogram,  as  the  dominat* 
period  of  the  maximum  phase  reduces  markedly  with  addition  of 
h-gner  terms.  Otherwise  the  pure  double-couple  (N  = 2)  seis- 
mogram includes  most  of  the  major  features.  Further  compu- 
tations verify  that  convergence  is  essentially  attained  at 
N - 6 • 

strat  lf0U“h  reC°rd  is  inoludeti  ^ figure  4.7  to  demon- 
trate  that  the  upper  mantle  response  used  in  the  computation 

1 not  frequency-independent  but  spreads  the  wave  form  some- 

(GS) *was°r  t^1S  reC°rd  3 constant  geometric  spreading  factor 

USS  t0  represent  the  "PPer  rather  than  the 

response  computed  with  a generalized  ray  theory  program. 

An  overall  view  of  the  effect  of  the  higher  order 
multipoles  (representing  propagation)  is  given  in  Figure  4.8 
Here  the  maximum  amplitudes  for  the  records  of  Figure  4 6 are 
normalized  to  the  value  for  4 = 45"  and  plotted  versus  azi- 
muth. For  a double-couple  the  amplitude  ratios  are  sin  24  as 
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Figure  4.7.  Theoretical  seismograms  at  several  levels  of 
truncation  (N)  of  the  multipolar  expansion. 
Otherwise  the  calculation  is  identical  to  that 
for  Figure  4.5a.  The  fourth  record  differs 
from  the  third  in  that  a constant  spreading 
factor  (GS)  was  used  to  represent  the  uppe~ 
mantle. 
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indicated  by  the  broken  line  on  the  figure.  (Note  that  the 
amplitude  for  the  N = 2 record  of  Figure  4.8  is  1.3  on  the 
scale  of  this  plot.)  If  mfa  is  plotted  in  this  way  the  con- 
trast with  a double-couple  pattern  is  somewhat  greater  due 
to  the  azimuthal  variation  of  T.  The  figure  shows  that  the 
higher  order  effects  have  a rather  small  influence  on  the 
average  amplitude  (or  n^)  for  a uniform  azimuthal  sampling 
but  cause  the  pattern  to  be  shifted  strongly  in  the  direc- 
tion of  fault  propagation. 

The  records  of  Figures  4.5  - 4.7  are  quite  simple; 
all  the  important  energy  is  in  the  P,  pP  and  sP  phases.  The 
lag  time  between  these  phases  is  determined  by  the  crustal 
structure,  which  otherwise  has  little  influence  on  the  seis- 
mogram. This  is  verified  in  Figure  4.9  where  seismograms 
are  shown  that  are  identical  to  those  in  Figure  4.5a  except 
that  the  source  region  crustal  structure  is  now  essentially 
a half space.  (The  properties  of  the  first  three  layers  of 

Table  4.1  were  changed  to  a,  6,  p = 7.9,  4.6,  3.3.)  Com- 
paring Figures  4.5a  and  4.9,  we  see  little  alteration  due  to 
reflections  from  velocity  contrasts  within  the  structure. 

The  pulse  distortion  due  to  the  layering  is,  of  course, 
strongly  dependent  or  frequency  but  the  long  periods  that 
dominate  these  records  do  not  "see"  velocity  contrasts  on 
the  scale  of  those  in  the  crustal  structure  used. 

The  amplitudes  on  the  seismograms  of  Figure  4.9  are 
relatively  constant;  variations  are  primarily  due  to  inter- 
ference between  the  direct  and  free  surface  reflected  phases. 
For  the  same  (in  terms  of  multipole  coefficients)  source  in 
a layered  crustal  model,  the  amplitudes  are  considerably 
more  varied  as  we  saw  in  Figure  4.5.  How  much  of  this  ampli- 
tude variation  is  dependent  on  the  source  parameter  changes 

(Table  4.3)  and  how  much  is  due  to  other  properties  of  the 
layered  model? 


L 
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Figure  4.8.  The  maximum  peak-to-peak  amplitude  versus  azi- 
muth for  the  strike-slip  source  at  25  km  focal 
depth.  The  amplitudes  are  from  the  records  of 
Figure  4.6  and  are  normalized  to  that  for 
0 = 45°.  The  closed  circles  represent  an  up- 
ward (compressional ) first  motion  and  the  open 
circles  are  for  downward  first  motion.  The 
broken  line  indicates  the  normalized  ampli- 
tudes expected  for  a simple  double-couple 
source.  The  circles  are  for  amplitude  ratios 
of  0.5  and  1.0. 
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Seismograms  like  those  of  Figure  4.5a  except 
that  the  first  three  layers  of  the  source 
region  crustal  structure  have  been  made  iden- 
tical to  the  fourth.  The  source  at  all  depth 
is  that  specified  by  (4.6). 
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Assume  that  we  can  measure  an  amplitude  on  the  seis- 
mogram that  is  proportional  only  to  the  direct  P wave  at  a 
fixed  frequency.  This  is  nearly  the  case  for  the  amplitude 
of  the  first  cycle  on  the  record  which  we  denote  by  b. 

Then 


b “ £ * ^ , 

where 

P — The  transfer  function  applied  to  the  displace- 
ments exiting  the  base  of  the  crust.  This  is 
a constant  for  all  seismograms  shown  here. 

T — The  transfer  function  representing  the  crust 
at  the  source. 

- The  source  generated  displacement  field  at  the 
period  at  which  b is  measured  and  along  the  ray 
that  propagates  to  the  receiver. 

From  (4.1)  and  (4.2)  the  U may  be  written  in  the 
simplified  form: 


(4.9) 


U * od  f(T,<j))  Aa(u) 
s s 


(4.10) 


where  ag  is  the  P wave  velocity  of  the  source  material, 
f(T,<j>)  is  the  radiation  pattern  effect  and  Aa(u>)  repre- 
sents the  multipole  coefficients.  For  all  the  seismograms 
of  Figures  4.5  - 4.7  and  4 . 9 we  have  fixed  v,  L/V  , v /8 
and  the  stress  drop  c(0).  The  A°(u>)  is  then  proportional 
to  strain  drop  for  which  the  variation  with  depth  is  indi- 
cated  in  Table  4.3.  The  superscript  o indicates  that 
0 The  f(i,<j>)  is  determined  from  radiation 

patterns,  like  those  of  Figure  4.3,  at  the  dominant  frequency 
for  the  b cycle.  The  takeoff  angle  for  the  f(T,<J>)  is  not 
constant  with  depth  but  varies  according  to  Snell's  law. 
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That  is. 


T “ sin"‘  as/c  , ( 

the  c being  fixed  by  the  upper  mantle  response. 

Bache,  et  aL  (1975)  showed  that  the  crustal  transfer 
function,  T,  reduces  to  approximately  a|/a*  where  a 
is  the  P wave  velocity  of  the  underlying  "basement"  mater- 
ial that  is  common  to  all  events.  Combining  all  these 
factors  we  see  that: 


b * ac  Aa(w)  . 

ft 


(4 


The  success  of  (4.12)  in  predicting  the  b amplitude 
for  the  seismograms  of  Figure  4.5  is  summarized  in  Table  4.4. 
The  A (W),  f(x ,<P)  and  the  predicted  and  measured  b are  all 
normalized  to  the  values  for  the  deepest  source.  The  T 
is  the  apparent  period  of  the  b phase  from  the  records  and 
is  seen  to  be  nearly  constant.  Radiation  patterns  at 
Tn  - 1.67  for  <p  = 45°  and  Tn  = 2.94  for  <p  = 135°  were  used 
to  determine  the  f(x,<f>).  This  value  was  then  adjusted  to 
the  apparent  period  of  the  measurement  by  multiplying  by 
(Tb/Tn)2'  thls  correction  being  motivated  by  the  slope  of 
the  spectrum  in  this  period  range  (see  Figure  4.4). 


Comparison  of  the  predicted  and  actual  normalized  b 
amplitudes  in  Table  4.4  indicates  that  (4.12)  gives  a reason- 
ably good  prediction  of  the  relative  amplitudes  on  the  body 
W?ve  seismograms.  The  differences  can  be  atrributed  to 
measurement  errors  and  the  fact  that  the  approximations  em- 
pxoyed  are  not  exact. 


The  scaling  given  by  (4.12)  should  be  taken  into  ac- 
count when  evaluating  the  dependence  of  teleseismic  ampli- 
tudes on  the  source  parameters.  We  have  employed  the  scalin 


• 11) 


• 12) 
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to  events  with  constant  v,  l/Vr,  Vr/«  and  o(0)  where  the 
multipole  coefficient  factor  A°(u)  , It  has  De£n  sug. 

gested  (e.g.,  Kanamori  and  Anderson,  1975)  that  strain  droo 
is  relatively  constant  for  large  events.  For  constant  strain 
drop  events  the  A <w>  is  replaced  by  AS(w>  which  is  in- 
variant for  the  events  compared  in  Table*4.4.  However,  in 
this  case  the  relative  amplitudes  on  the  seismograms  are 

scaled  by  ratios  of  ug  and  the  normalized  amplitudes  in 
the  table  remain  the  same. 

Table  4.4.  Scaling^of  b Amplitudes  for  the  Records  of 


Depth 

(km) 

as 

(km/sec) 

Aa(a<°)  . 
100  bars) 

f (6,0) 

Predicted 

b 

Actual 

b 

T_b 

o 

in 

■tr 

II 

-©- 

10 

6.0 

2.03 

0.53 

0.35 

0,37 

1.63 

15 

6.0 

2.03 

0.61 

0.41 

0.48 

1.75 

25 

6.7 

1.43 

0.64 

0.47 

0.55 

1.67 

35 

7.9 

1.00 

0.98 

0.98 

0.90 

1.72 

45 

7.9 

1.00 

1.00 

1.00 

1.00 

1.74 

<f>  = 135° 

10 

6.0 

2.03 

0.64 

0.43 

0.45 

2.90 

15 

6.0 

2.03 

0.55 

0.37 

0.38 

2.70 

25 

6.7 

1.43 

0.69 

0.51 

0.55 

2.80 

35 

7.9 

1.  00 

1.05 

1.05 

1.08 

3.07 

45 

7.9 

1.00 

1.00 

1.00 

1.00 

3.00 
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4.6  CONCLUDING  REMARKS 

If  we  restrict  attention  to  body  wave  recordings  at 
distances  beyond  the  upper  mantle  triplication  and  from 
stations  where  the  local  crustal  reverberations  are  small, 
most  of  the  travel  path  effects  are  minimized  and  the  record 
complexity  can  be  attributed  to  the  complexity  of  the  source. 
This  observation  is  born  out  by  at  least  some  teleseismic 
recordings  of  large  underground  explosions.  However,  even 
when  the  travel  path  effects  are  minimized,  most  earthquake 
recordings  remain  quite  complex.  In  this  sense  the  syn- 
thetic seismograms  presented  here  are  too  simple  to  be 
realistic.  An  actual  earthquake,  especially  one  with  a 
fault  length  near  10  km,  can  hardly  be  a single  event  char- 
acterized by  a single  stress  drop,  rupture  velocity,  etc. 

In  fact,  we  would  expect  such  an  earthquake  to  be  composed 
of  many  smaller  events,  each  with  its  own  fault  parameters. 
From  this  point  of  view  our  single  source  model  is  too 
idealistic  to  provide  a detailed  synthesis  of  an  actual  earth- 
quake recording.  But  this  was  not  really  our  objective. 

The  development  of  a capability  for  simulating  the 
earthquake  source  and  propagation  of  the  resulting  seismic 
waves  is  a step-by-step  process  in  which  increasing  levels 
of  complexity  are  added  as  the  need  arises.  Techniques  for 
including  more  detail  in  both  the  source  and  propagation 
models  are  being  considerably  improved  and  the  computational 
methods  used  in  this  paper  were  designed  to  link  the  two. 

At  least  one  important  "higher  order"  source  effect  is  con- 
veniently included  in  the  Archambeau/Minster  earthquake 
model  used  here;  that  is,  the  effect  of  finite  velocity 
rupture  propagation.  The  far-field  radiation  from  such  a 
source  is  quite  different  from  that  from  a pure  double- 
couple source  representation  as  we  saw  in  the  patterns  of 
Figure  4.3.  For  teleseismic  body  wave  observations,  the 
consequences  include  a pronounced  shift  of  the  teleseismic 
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body  wave  radiation  pattern  in  the  direction  of  rupture 
propagation  (Figure  4.0). 

A scaling  Jaw  for  teleseismic  body  wave  amplitudes, 
subject  to  certain  similarity  constraints  on  the  source,  is 
an  interesting  consequence  of  the  seismogram  computations 
presented.  The  first  constraint  is  that  Poisson's  ratio  be 
fixed  for  the  source  regions  being  compared.  Then  if  the 
ratio  of  fault  length  to  rupture  velocity  (L/VR)  and  rupture 
velocity  to  shear  velocity  (L/8)  are  fixed, 

a = a4  f(T,$)  A{(j)  / 

* ' 

where  a represents  the  amplitude  of  either  the  direct  P 
or  S wave.  The  f(f,<j>)  is  the  radiation  pattern  effect  which 
depends  on  the  azimuth  (<J>)  and  takeoff  angle  (x)  for  the  ray 
contributing  to  the  a amplitude.  The  A(w)  represents 
the  source  strength  (the  multipole  coefficients)  and  it  is 
to  this  term  that  the  similarity  constraints  apply.  The 
dependence  of  A(uj)  on  the  source  parameters  is  discussed 
by  Archambeau  and  Minster  (1976).  One  important  property 
-(that  A(w)  is  directly  proportional  to  strain  drop 
(S  ).  The  a4  factor  includes  a2  for  converting  A(ca) 
to  displacement  and  a2  which  accounts  for  the  propagation 
from  the  source  material  into  the  underlying  mantle  material 
that  is  assumed  to  be  fixed  for  the  events  being  compared. 
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V. 


THE  MULTIPLE  ARRIVAL  RECOGNITION  SYSTEM 


5* 1 INTRODUCTION 

of  th  Th,iS  SeCti°n  °f  the  rePOrt  iS  deVOted  t0  * description 
the  mode  of  operation  of  the  Multiple  Arrival  Recognition 

ystem  (mars)  computer  code  and  the  application  ol  this  code 

explosion  seismograms.  The  MARS  code  is  proving  to  be 

particularly  effective  for  detecting,  isolating  and  timing 
the  various  seismic  phases  (P  , P*  p - 9 

G ' n' 


that  are  recorded  on  event  seismogram^ 'in ne«'r!gIoLl 
and  reqional  y onai 


, , ah  uie  near-reqional 

and  regional  distance  ranges.  This  i 

is  finrUn-,  g h Signal  analysis  capability 

is  finding  immediate  application  <-o  . 

application  to  the  problem  of  depth  de- 

terminations  for  shallow  focus  events.  More  accurate  event 
epth  determinations  can  be  obtained  by  incorporating  P and 

and  bv  ^riVal  tlmES  in  sta"dard  hypocentral  location  routines 
y he  combined  use  of  spectra  of  the  isolated  body  waves 
or  identification  of  spectral  nulls  caused  by  the  inter 
erence  of  direct  and  free-surface  reflected  phases.  Future 
t a S on  the  MARS  code  should  involve  a systematic  analysis 
Of  three-component  data  recorded  at  near-regional,  regional 
E eseismlc  distances  from  Eurasian  events. 

As  reported  in  previous  reports  written  under  this 

relat^l  ~ ' 1975;  Savino,  et  ah,  1975),  the 

the  variable  f °tiVe  referred  to  as 

em a "a9nitude  (VFM)  discriminant,  is 

embodied  in  the  MARS  code,  m chapter  VI  of  this 

T Will.deSCribe  the  of  the  application  of 

served  lmi"ant  t0  ^"thetic  seismograms  and  to  ob- 

n 112S:Z7\T  eamqUa,te  and  eXPl0Si°"  Populations 
urasia  and  North  America. 


5.2 


DESCRIPTION  of  MARS 


is  the  use  Ttral  featUre  °£  the  "MS  3^nal  analysis  code 
use  of  narrow  band  frequency  filters  to  break  up  or 
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decompose  a time  series  consisting  of  signal  plus  noise  or 
just  noise  into  a set  of  quasi-harmonic  modulated  "signals". 
This  set  of  filtered  signals,  one  for  each  filter  center 
frequency,  can  then  be  used  to  determine  the  energy  arrival 
time  (the  group  time,  tg)  and  amplitude  of  the  signal  for 
each  filter  frequency  by  analysis  of  the  time  modulation  of 
the  filter  outputs.  Further,  both  the  instantaneous;  phase 
and  frequency  of  the  individual  filter  outputs,  that  is  the 
apparent  phase  and  frequency  of  the  quasi-harmonic  filter  out- 
put as  a function  of  time,  can  be  determined  quite  simply. 
Thus,  if  two  components  of  the  same  vector  wave  field  are 
analyzed  in  this  fashion,  the  difference  in  the  instantaneous 
phase  between  the  two  components  is  a direct  measure  of  the 
polarization  of  the  wave  field. 

Thus,  the  decomposition  of  the  signal  wave  train,  pos- 
sibly composed  of  many  individual  signal  pulses,  into  quasi- 
harmonic signals  provides  the  means  of  determining  arrival 
time,  amplitude  and  wave  polarization  all  as  functions  of 
frequency.  This  then  is  the  basic  signal  information  that 
can  be  used  to  detect  a given  type  of  signal  in  terms  of  its 
dispersion  and  polarization  characteristics  and  to  obtain 
its  spectrum  as  well  as  its  time  and  amplitude  relation  with 
respect  to  other  signals  present  in  a complex  wave  train. 

Figures  5.1  and  5.2  are  two  different  versions  of  the 
flow  of  operations  in  the  MARS  code.  Figure  5.1  gives  a 
verbal  outline  while  Figure  5.2  summarizes  the  key  mathemati- 
cal operations  performed  in  this  signal  analysis  code.  More 
detailed  developments  of  the  theory  and  operation  of  MARS 
were  presented  in  Bache,  et  al^  (1975)  and  Savino,  et  al^ 

(1975) . ^ ” 

Seismic  data  are  read  into  MARS  in  the  form  of  a time 
series  (Figures  5.1  and  5.2)  generally  of  about  500  - 2000 
points  in  length.  The  data  are  then  optionally  detrended. 
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Plot  of  Envelopes  and  Sroup  Arrival  T'oa 


Capita  Polarization  r 11 tar  Functioni  for  Body 
<P>S)  and  Surfaca  CLR.LQJ  Vkvaa;  for  Exanpla, 
for  P wave*  - 

p3’  - *•"  ["S’]  «■"  tan'i(^r)  ■ a 

vtwxrt  n,  m art  Xntagars,  2 an  Bnar^anca  Angla. 

(k) 


Flowchart  indicating  principal  mathematical  opera- 
tions embodied  in  th  ''a.s  code. 
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mean  removed  and  tapered  at  the  tail  end  by  a cosine  bell. 
MARS  then  selects  the  smallest  power  of  two  which  is  greater 
than  the  number  of  points  input  and  performs  a discrete 
Fourier  transform  using  the  algorithm  c£  Cooley  and  Tukey 

(1965).  Both  the  original  time  serxes  and  the  spectrum  are 
plotted  for  examination. 


Referring  to  Figures  5.1  and  5.2,  the  signal  is  next 
filtered  in  the  frequency  domain  by  multiplication  with  a 
narrow-band  cusp-shaped  filter.  This  particular  filter  form 
is  selected  to  satisfy  two  goals:  (1)  minimum  width  in  the 

frequency  domain,  and  (2)  maximum  ripple  suppression  in  the 
time  domain.  A well-known  consequence  of  the  Uncertainty 
Principle  (Sampling  Theorem  to  electrical  engineers)  is  that 
one  cannot  simultaneously  satisfy  these  two  goals  to  arbitrary 
precision.  The  filter  employed  in  MARS  was  selected  for  its 
optimal  time  and  frequency  domain  characteristics  within  the 
limits  of  the  Uncertainty  Principle. 


Once  the  signal  has  been  narrow-band  filtered,  MARS 
corrects  it  for  the  appropriate  instrument  response;  the 
filtered  signal  transform  is  divided  by  the  instrument  trans- 
fer function.  The  resulting  complex  spectrum  is  then  inverse 
Fourier  transformed  into  the  time  domain,  to  produce  what  will 
hereafter  be  referred  to  as  the  filtered  signal. 


The  narrow-band  filtered  sig  a wi  1 appear  as  a quasi- 
sinusoidal  carrier  wave  contained  wit  .in  a smooth  envelope. 
The  next  step  in  the  MARS  code  (Figures  5.1  and  5.2)  is  to 
construct  the  envelope  function  by  means  of  the  Hilbert  trans- 
form. This  method  is  followed  in  MARS:  the  transform  of  the 

filtered  signal  is  multipled  by  -i  sgn  («,)  and  then  brought 
to  the  time  domain  by  an  inverse  transformation.  The  maximum 
of  the  envelope  function  is  utilized  in  the  VFM  discriminant 
W He  the  instantaneous  frequency  and  phase  are  stored  for  sub 

sequent  use  in  polarization  filtering  with  additional  componen 
of  ground  motion. 
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The  narrow-band  filtering  procedure  can  be  performed 
on  a particular  component  seismogram  (time  series)  at  a num- 
ber of  different  frequencies  within  some  band  of  interest. 
Correlation  of  the  resulting  envelope  functions  indicates 
the  arrival  times  of  the  various  frequency  components.  Then, 
for  example,  the  group  velocity  dispersion  can  be  deduced 
from  the  variation  of  the  wave  group  arrival  time  (t  ) with 
filter  center  frequency  (fj . An  application  of  mar!  for 
determining  Rayleigh  wave  dispersion  along  a North  American 
travel  path  was  discussed  by  Barker,  et  al.  (1976) . 

In  the  event  that  multicomponent  ground  motion  data 
are  available,  polarization  filters,  designed  to  pass  ground 
motion  having  the  appropriate  particle  motion,  can  be  applied 
to  the  filtered  data  for  the  purpose  of  identifying  and 
separating  waves  of  differing  polarization.  Polarization  fil- 
tering in  MARS  is  always  performed  on  time  series  which  have 
previously  been  narrow-band  filtered.  The  rationale  for  this 
procedure  is  that,  in  general,  the  spectral  content  of  the 
signal  (e.g.,  body  wave)  differs  from  that  of  the  recorded 
background  noise.  As  a result,  the  application  of  filters  of 
appropriate  frequencies  can  usually  enhance  the  signal-to- 
noise  ratio  thereby  partially  isolating  frequency  components 
of  the  signal  for  further  processing. 


~ • 3 mars  ANALYSIS  OF  THE  NTS  EXPLOSION  CHARTREUSE 

On  May  6,  1966,  at  15:00:00.1  GMT,  the  underground 
explosion  code-named  Chartreuse  was  detonated  in  the  Pahute 
Mesa  area  at  the  Nevada  Test  Site  (NTS) . Digitally  recorded 
seismograms  for  Chartreuse  were  obtained  from  two  stations 
of  the  Long-Range  Seismic  Measurements  (LRSM)  network  for 
processing  with  the  MARS  code.  These  stations  and  epicentral 
distances  are:  Mina,  Nevada  (MNNV)  at  201  km  and  Kanab,  Utah 

(KNUT)  at  312  km. 
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We  are  primarily  interested  in  detecting  and  identify- 
ing the  body-wave  phases  that  are  not  commonly  noted  in  near 
regional  explosion  seismograms.  To  detect  and  identify  shear 
vertical  phases  hidden  in  dominating  Rayleigh  wave  signals 
requires  the  utilization  of  polarization  filtering. 

The  instantaneous  phase  <f>n  determined  from  variable 
frequency  narrow-band  filtering  is  input  directly  to  the 
polarization  filter  routine  of  the  MARS.  The  linear  polariza- 
tion fraction  can  be  defined  using  the  instantaneous  phase 
difference  between  the  vertical  (z)  and  radial  (r)  seismic 

components.  Thus,  the  linear  polarization  fraction  (P  ) is 
defined  as  ** 

PL(t)  = cos  [0z(t)  - 4>r(t)]  , 

and  the  circular  polarization  fraction  as 
Pc(t)  = sin  [<J>z  (t)  - <Dr(t)]. 

Therefore,  it  is  clear  that 
PL  + Pc  " 

and  the  two- are,  in  a sense,  mutually  exclusive  (Savino . et 
al. , 1975).  ‘ — 

For  z-r  data,  different  types  of  wave  motion  can  be 
identified: 

P-wave:  Pj.  > 0 and  PL  > |p  | . 

SV-wave:  PL  < 0 and  Jp  | > |p  | ; 

Xj  C 

Rayleigh  wave:  Pc  < 0 and  | ?c | > |p  | . 

These  types  of  motion  can  be  displayed  in  time  by  multiplying 
the  polarization  fraction  having  the  appropriate  constraints 
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by  the  quasi-harmonic  modulated  waves  resulting  from  the 
narrow-band  filtering.  Such  results  can  then  be  displayed 
in  the  time  verrus  frequency  plane  to  provide  arrival  time 
estimations  by  alignment  of  non-dispersed  P and  S waves  and 
indicated  dispersive  qualities  of  the  Rayleigh  waves.  Re- 
finement of  arrival  time  estimations  is  achieved  by  forming 
the  envelope  functions  (S)  and  averaging  the  envelope  peaks 
over  many  frequencies.  The  envelope  functions  are  formed  in 
the  following  manner. 

For  the  P wave 

Sp(t)  = max  (PL(t),0)  • (t)  + S*  (t)  . 

For  the  SV-wave 

Sgv(t)  = max  (-PL(t),0)  • yjs2z(t)  + S*  (t)  , 
and  for  the  Rayleigh  wave 

S (t)  = max  (-P  (t),0)  • ^ S2  (t)  + S2  (t)  , 
c z r 

where  Sz  and  are  narrow-band  filtered  time  series  of  the 
vertical  and  radial  seismic  components. 

Figure  5.3a  shows  the  recorded  vertical  ground  displace- 
ment at  Kanab,  Utah  (KNUT) . As  a consequence  of  the  short 
distance  (312  km)  between  source  and  receiver,  the  record  is 
very  complicated,  with  many  arrivals  compressed  into  a fairly 

time  interval.  The  time  t = 0 on  this  record  is  actually 
43.1  seconds  after  detonation  time. 

Simple  visual  inspection  reveals  that  for  the  first  20- 
25  seconds,  the  record  is  dominated  by  motion  at  a relatively 
high  frequency,  the  amplitude  of  which  then  decreases  as  the 
record  is  taken  over  by  lower- frequency  oscillations.  In 
fact,  at  approximately  27  seconds  into  the  record,  two  distinct 


frequency  components  can  be  clearly  seen.  It  would  appear, 
from  observation,  that  some  improvement  in  phase  separation 
could  be  made  by  narrow-band  and  polarization  filtering. 

The  vertical  and  radial  components  of  the  Chartreuse- 
KNUT  record  were  processed  by  the  MARS  code,  with  an  eye  to 
phase  separation.  ne  rapidly  varying  oscillations  at  early 
times  appear  to  be  at  approximately  1.5  Hz.  When  this  fre- 
quency is  used  for  the  narrow-band  filter,  the  resulting  P- 
and  SV-waves  are  as  shown  in  Figures  5.3b  and  5.3c,  respectively. 

Formulae  for  the  travel  times  of  certain  body  wave 
phases  were  given  by  Gutenberg  (1944)  and  modified  by  Evernden 
(1967)  who  shows  that  for  the  path  from  NTS  to  KNUT  a Moho  P 
velocity  of  7.8  to  7.9  is  appropriate.  For  the  travel  path 
of  interest,  the  travel  times  for  P-waves  are  as  given 
below: 


t(pn) 

- t(0) 

= X + 

= 46.2  sec  , 

t(p*> 

- t (0) 

= 4.4 

+ (T794  = 49-4  sec  ' 

t(p‘) 

- t(0) 

= 1.2 

+ o?  " 52-8  sec  ' 

t<p») 

= t(0) 

■ = 55.9  sec  . 

5.58 

In  these  equations,  X = 6.2  is  the  surface-to-Moho  travel  time, 
and  A = 312  km.  The  four  computed  times  are  indicated  in 
Figure  5.3b;  it  must  be  kept  in  mind  that  the  origin  is  shifted 
by  43.1  seconds.  It  is  clear  that  MARS  has  identified  P wave 
phases  with  arrival  times  approximately  the  same  as  in  the 
tables.  The  remaining  bursts  of  1.5  Hz  P— wave  energy  are  not 
currently  identifiable,  and  are  presumably  due  to  more  compli- 
cated travel  times,  involving  multiple  reflections. 
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Gutenberg  (1944)  gives  travel  time  equations  for  cor- 
responding shear  body  phases.  These  times  work  out  to  be: 


t(S  ) 
n 

- t (0) 

= Y + 

4.31  " 80*9  sec 

t (S*) 

- t(0) 

= 6.9 

+ 4>1Q  = 83.0  sec 

t(S‘) 

- t(0) 

= 2.1 

+ 3.65  " 87,6  sec 

t(S2) 

- t(0) 

A 

" 1728 

■ = 95.1  sec 

where  Y 8.5  seconds  and  A = 312  km.  The  appropriate  times 
are  indicated  in  Figure  5.3c,  where  they  clearly  correspond 
to  recognizable  S-wave  arrivals.  Additional  wavepacket.s  are 
visible,  some  at  earlier  times,  but  these  are  presumably 
attributable  to  P-S  converted  waves  which  travel  most  of 
the  distance  at  the  longitudinal  wave  speed. 

The  low  frequency  oscillations  which  are  visible  in 
Figure  5.3a  can  be  brought  out  by  the  application  of  an  0.9 
Hz  narrow-band  filter.  When  this  is  done,  the  wave  appears 
quite  distinctly,  the  F*  disappears,  while  the  P1  and  P2  merge 
into  one  wavepacket.  The  0.9  Hz  P-waves  are  not  illustrated 
' since  they  do  not  show  anything  not  previously  seen. 

For  0.9  Hz,  the  projected  Rayleigh  waves  are  as  shown 
in  Figure  5.3d;  several  clear  arrivals  can  be  seen  between 
record  times  of  25  and  55  seconds.  Given  the  distance  of 
312  km,  and  uhe  t = 0 offset  of  43  seconds,  these  wavepackets 
have  been  traveling  at  between  3.2  and  4.6  km/sec.  These 
group  velocities  are  consistent  wrth  the  values  for  higher 
mode  Rayleigh  waves  observed  by  others  for  the  Basin  and 
Range  region,  Alexander  (1963). 

Figure  5.4a  shows  the  vertical  ground  displacement  at 
Mina,  Nevada  (MNNV) . In  general  this  record  is  similar  to 
that  shown  for  KNUT  (Figure  5.3a).  However,  this  station  is 
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located  201  kilometers  from  the  Chartreuse  epicenter,  which 
is  considerably  closer  to  the  explosion  than  KNUT  (312  km) . 
Thus,  similar  individual  seismic  signals  identified  at  the 
KNUT  station  will  arrive  earlier  and  closer  together  in  time. 
Consequently,  separation  of  these  individual  seismic  signals 
can  only  be  accomplished  with  the  use  of  relatively  higher 
frequency  narrow-band  filters. 

In  the  following  figures,  5.4b,  5.4c  and  5.4d,  we  show 
the  polarization  filter  outputs  for  the  frequency  (2.8  Hz 
body  waves  and  0.9  Hz  Rayleigh  waves)  that  best  delineate  in 
time  the  individual  seismic  signals.  At  this  station  we  are 
unable  to  identify  a P*  or  S*  (Figures  5.4b  and  5.4c). 
Gutenberg  (1944)  and  Herrin  (1961)  travel-time  tables  indicate 
that  for  this  distance  P*  should  arrive  about  one  second  later 
than  the  Pn-  On  the  other  hand,  Herrin  (1968)  travel-time 
curves  indicate  that  P*  should  arrive  several  tenths  of  a 
second  earlier  than  Pn«  Further,  for  the  total  series  of 
narrow-band  filters  applied,  we  observed  no  clear  seismic 
P arrival  between  the  Pn  and  P^  seismic  phase. 

As  for  KNUT,  we  show  the  projected  Rayleigh  waves  at 
0.9  Hz  for  MNNV , and  as  before  the  croup  velocities  of  the 
major  wave  packets  are  consistent  with  the  values  for  higher. 

mode  Rayleigh  waves  ouserved  by  others  for  the  Basin  and  Range 
region. 

Rayleigh  wave  motion  dominates  the  later  (in  time) 
portion  of  both  the  MNNV  and  KNUT  seismograms;  however,  polari- 
zation filtering  clearly  defines  the  important  SV  phases  which 
would  be  most  difficult  if  not  impossible  to  pick  from  the 
original  seismograms  (Figures  5.3a  and  5.4a). 

Comparison  of  the  predicted  P and  S travel-times  to 
the  corresponding  observed  travel- times  are  given  in  Table 

5.x.  Differences  between  predicted  and  observed  times  range 
up  to  0.7  second. 
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Table  5.1.  Predicted  Versus  Observed  Travel-Times 


Travel-Times  f nr 

KNUT 

Phase 

Predicted 

(sec) 

Observed 

(sec) 

Residual  (sec) 
(Predicted-Observ^  ^ 

P 

n 

46.2 

46.4 

• 

-0.2 

P* 

49.4 

48.7 

+0.7 

P1 

g 

52.8 

52.8 

0.0 

P2 

g 

55.9 

56.2 

-0.3 

s 

n 

80.9 

80.4 

+0.5 

S* 

83.0 

82.7 

+0.3 

sg 

87.6 

88.1 

-0.5 

95.1 

95.2 

-0.1 

Travel-Times  f nr 

MNNV 

P 

n 

32.0 

32.3 

-0.3 

P1 

g 

34.4 

34.1 

+0.3 

P2 

g 

36.0 

36.4 

-0.4 

Sn 

55.1 

55.4 

-0.3 

sg 

57.2 

57.3 

-0.1 

s2 

g 

61.3 

61.0 

+0.3 
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Thus,  with  the  use  of  MARS  we  have  detected,  identi- 
fied, isolated  and  accurately  times  P,  s and  Rayleigh  seismic 
phases  at  a small  epicentral  distance.  Especially  noteworthy 

IS  the  detection  of  Sn  since  this  phase  is  rarely  observed 
from  explosions. 


discussion  of  depth  estimation 


5*4.1  Location  Programs 

If  very  accurate  estimates  of  the  origin  time  are 
available,  one  can  generally  determine  the  depth  of  focus 
with  considerably  less  error  than  when  neither  origin  times 
nor  focal  depth  are  constrained.  In  many  cases,  one  can 
use  P and  S travel-times  together  to  improve  estimates  of 
location  and  depth,  though  this  fact  has  not  been  exploited 
in  widely  used  location  programs. 


One  method  of  determining  the  origin  time  of  a seismic 
event  is  to  measure  the  arrival  times  of  P and  S waves  that 
follow  the  same  propagation  path  from  source  receiver  (e.g., 
t(Pn),  t(Sn);  t(pg),  t(Sg);  and  etc.).  Several  stations  at 
varying  distances  are  required  such  that  when  plotting  P 
versus  S arrival  times,  the  point  t(S)  - t(P)  = 0 is  well 
determined.  This  time  will  be  the  event  origin  time. 


Evernden  (1969)  performs  an  exercise  utilizing  P and 
Sg  arrival  times  for  the  Bilby  and  Aardvark  NTS  explosions 
and  determines  origin  times  within  0.5  seconds  of  the  real 
times.  He  also  points  out  that  the  errors  in  picking  the 
pg  and  sg  times  are  commonly  on  the  order  of  2 seconds.  Ap- 
plication of  MARS  for  determination  of  various  P and  S arrival 
times  will  significantly  reduce  the  error  and  the  number  of 
stations  required.  This  allows  a more  accurate  estimate  of 
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or^9in  time.  Evernden  (1969)  indicates  that  if  the  origin 
time  is  constrained  with  an  error  of  0.5  seconds,  the  loca- 
tion computation  gives  depth  estimates  of  0 +4.5  km  for 
NTS  explosions.  The  travel-time  curves  used  by  Evernden  in 
the  location  computation  were  those  of  Herrin  (1961) . These 
are  not  entirely  accurate  for  the  western  United  States, 
but  even  if  they  were,  the  depth  errors  would  have  been  only 
slightly  less  so.  Very  accurate  travel— time  curves  are 
required  for  accurate  depth  estimations  when  using  uncon- 
s**ra^ne<^  origin  times  and  depths  in  the  location  computation. 

Figure  5.5  is  a plot  of  the  P and  S arrival  times  de- 
termined and  discussed  previously  for  MNNV  and  KNUT.  A mean 
line  determined  by  a least  squares  fitting  to  the  origin 
time  of  14:59:59.8  or  0.3  seconds  earlier  than  the  actual 
origin  time.  The  standard  deviation  of  the  points  normal 
to  the  line  is  + 0.3  seconds.  These  results  are  remarkable. 
However,  we  must  keep  in  mind  that  only  two  stations  were 
used  and  the  excellent  results  may  be  fortuitous.  More 
station  data  at  differing  distances  would  provide  added  con- 
fidence in  the  origin  time  estimate. 

In  conclusion,  sv cstantially  improved  depth  estimates 
can  be  obtained  by  using  MARS  to  pick  accurate  arrival  times 
for  P and  S phases  from  seismograph  recordings  at  distances 
of  150  to  1500  km.  With  the  accurate  arrival  times  we  can 
closely  constrain  the  event  origin  time  which,  in  turn,  is 
used  to  constrain  the  depth  estimate. 

5-4.2  Spectra  for  Narrow-Band  Filtering 

When  the  location  computation  indicates  that  the 
depth  is  shallow  (0  + 4 or  5 km) , the  spectra  of  the  various 
body  waves  that  are  identified  and  isolated  with  the  use  of 
the  MARS  can  be  utilized  to  refine  the  depth  estimation. 
Spectral  nulls  will  occur  due  to  the  interference  in  the  time 


Arrival  Times  of  S(t  ) 


Figure  5.5. 
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series  of  the  depth  phase  (free-surface  reflection)  with  the 
more  direct  downgoing  wave. 

Phase  interference  of  this  type  may  be  represented  by 
the  following: 

z (t)  = y (t)  - By  (t  + t) 


where 

t = delay  time  associated  with  the  surface  reflection, 
y(t)  * initial  source  time  function, 

8 = the  surface  reflection  coefficient  (6  <_  1.0). 

The  power  spectrum  of  z (t)  is  then: 


P(u>) 


I / 


z (t)  e 


-loot 


dw 


which  becomes 

P(w)  = |Y(w)|2  [1  + B2  - 28cosojt]  . 

Thus,  nulls  in  the  spectrum  will  occur  1/x  Hz  apart. 

In  Figure  5.6,  we  show  the  amplitude  spectra  for  the 

seismic  phases  P^,  P* , and  P^  as  obtained  from  the  series  of 

narrow-band  filters  applied  to  the  Chartreuse-KNUT  record. 

We  would  expect  that  the  delay  times  t would  increase  with 

the  increase  of  distance  that  the  reflected  wave  must  travel 

due  to  the  increasing  angle  of  incidence  corresponding  to 

each  of  these  phases  (P  . P*  and  P2).  This  should  be  es- 

n g 

pecially  true  for  P2  as  compared  to  P since  P2  phase  travels 

9 n 9 

in  the  upper  crust  layers.  In  Figure  5.6  we  observe  that 
x does  increase  as  expected.  This  type  of  comparative  spectral 
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analysis  greatly  increases  the  confidence  that  can  be  placed 

on  the  selection  of  the  spectral  nulls  for  depth  estimations 

of  shallow  events  observed  at  one  station.  The  inclusion  of 

other  station  data  will  have  the  same  effect.  In  fact,  from 

stations  having  common  seismic  phase  spectra,  summation  of 

these  spectra  will  enhance  the  spectral  nulls  due  to  depth- 

phase  interference  and  tend  to  cancel  varying  path  and  re- 
ceiver  effects . 


5 * 5 SUMMARY  AND  RECOMMENDATIONS 

We  have  demonstrated  in  part  the  utility  of  the  MARS 
code  from  the  seismic  analyst's  point  of  view.  The  specific 
problem  of  seismic  event  depth  estimation  was  addressed.  In 
this  respect  the  MARS  capability  to  detect,  isolate,  identify 
and  accurately  measure  travel  times  and  amplitudes  of  seismic 
fhases  making  up  an  extremely  complicated  record  can  be  used 

m a simple  and  effective  way  to  obtain  accurate  estimates 
Of  event  depth  of  focus. 

In  particular,  the  use  of  accurately  determined  arrival 
times  for  various  P and  S waves  that  follow  the  .-'me  propaga- 
tion path  from  source  to  receiver  to  obtain  event  origin 
times  will  yield  improved  depth  estimations  in  the  location 
computation.  For  those  very  shallow  events,  such  as  under- 
ground explosions,  the  use  of  various  individual  body  wave 
spectra  to  clearly  identify  spectral  nulls  should  virtually 
solve  the  problem  of  depth  determinations  for  shallow  seismic 

events,  especially  so  if  data  from  several  recording  sites 
are  combined. 

Future  work  should  be  concerned  with  a systematic  test 
of  the  discussed  procedures  on  a suite  of  Eurasian  events  re- 
corded at  stations  located  within  near  regional,  regional  and 
small  teleseismic  distances. 
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VI.  VARIABLE  FREQUENCY  MAGNITUDE  DISCRIMINANT 


6 . 1 INTRODUCTION 


An  important  aspect  of  the  research  program  supported 
by  the  subject  contract  has  been  the  development  and  testing 
of  optimal  procedures  for  discrimination  between  earthquakes 
and  both  single  and  multiple  underground  explosions.  As 
previously  reported  (Baohe,  et  aU . 1974;  Savino  and  Archambeau, 
1974;  Baohe,  et  al.,  1975;  Savino,  et  al..  1975),  this  program 
has  resulted  in  the  formulation  of  a very  effective  variable 
frequency  magnitude  (VFM)  discriminant  which  exploits  spectral 
differences  between  earthquakes  and  explosions. 


To  date  our  attention  has  focussed  on  the  application 
of  the  VFM  technique  to  short  period  body  waves  recorded 
from  events  at  teleseismic  distances.  A well-known  limita- 
tion to  the  applicability  of  the  much  heralded  M -„i  discrimi- 
nant results  from  the  inability  of  presently  existing  seismic 
networks  to  detect  surface  waves  from' low-yield  explosions. 

For  example,  Basham  and  Whitham  (1971)  estimated  that  the  de- 
tection threshold  for  20  second  Rayleigh  waves  from  explosions 
was  0.5  to  0.7  magnitude  units  higher  than  the  p wave  thres- 
hold for  the  same  event  population.  In  the  case  of  the  VFM 
discriminant,  obviating  the  requirement  for  recording  surface 
waves  from  small  explosions  results  in  an  extended  magnitude 
range  over  which  this  discriminant  can  be  applied.  Based 
on  the  North  American  and  Eurasian  event  populations  examined 
so  far,  complete  separation  of  earthquakes  and  explosions  is 
achievec.  with  the  VFM  technique  down  to  event  magnitudes,  m„, 
in  the  range  4.0  to  4.5.  b 


Following  a brief  description  of  the  physical  bap is 
for  the  VFM  discriminant,  the  computational  procedure  followed 
in  applying  this  technique  will  be  outlined.  Next  the  results 
of  the  application  of  the  VFM  technique  to  several  large 
populations  of  explosions  and  shallow  and  deep  focus 


earthquakes  recorded  at  the  LASA,  NORSAR  and  Yellowknife, 
Canada,  arrays  will  be  discussed.  Then  the  effectiveness 
of  the  VFM  technique  for  discriminating  a simulated  multiple 
explosion  scenario  will  be  tested.  This  section  of  the 
report  will  conclude  with  the  application  of  the  VFM  dis- 
criminant to  synthetic  seismograms  for  earthquakes  and 

explosions,  and  a comparison  of  the  synthetic  and  observed 
VFM  results. 


6-2  PHYSICAL  BASIS  FOR  THE  VFM  DISrRTMTMaMm 

The  variable  frequency  magnitude  discriminant,  originally 

proposed  by  Archambeau,  et  ah  (1974,,  exploits  spectral  dif 

erences  between  underground  explosions  and  earthquakes.  The 

approach  employed  makes  use  of  spectral  body  wave  magnitudes 

computed  from  narrow-band  filter  outputs  at  frequencies 

generally  in  the  range  0.4  Hz  to  5 h»  a.  , _ 

u. « hz  to  5 Hz,  the  actual  frequencies 

being  dependent  on  the  particular  source-to-receiver  propaga- 
tion path  studied.  These  spectral  magnitudes  are  denoted  as 
y f>-  m this  report  we  will  consider  body  wave  magnitudes 
at  different  frequencies  (i^  (f , , , ,^<f  f ) ) for  discrimination. 

Differences  in  the  far-field  short  period  P-„ave  spectra 
or  earthquakes  and  underground  explosions  have  been  reported 
by  many  investigators  (Lacoss,  1969;  Wyss,  et  , 1971;  Molnar, 
1971).  m general  the  observations  indicate  that  if  explo- 
sion and  earthquake  spectral  levels  are  comparable  at  fre- 
quencies lower  than  1 Hr  (e.g.,  0.5  Hr),  then  at  higher  fre- 
quencies (e.g.,  2 HZ)  the  explosion  spectrum  is  enriched 
relative  to  that  for  the  earthquake.  These  observations  led 
to  the  rormulation  of  several  different  discrimination  techni- 
ques designed  to  exploit  these  spectral  differences;  spectral 
ratio  (Lacoss,  1969),  and  higher  moments  of  frequency 
(Anglin,  1971;  Manchee,  1972). 

The  combined  effects  of  the  free  surface  reflection, 
smaller  source  dimension,  and  peaked  source  function  have 
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been  invoked  to  explain  the  shape  of  short-period  explosion 
spectra  relative  to  earthquake  spectra.  The  observed  degrada- 
tion of  the  lower  frequency  portion  of  explosion  spectra  has 
been  attributed  to  the  superposition  of  the  direct  P wave  and 
the  free  surface  reflected  wave,  pP  (lagged  in  time  and  of 
opposite  polarity  to  P) , which  results  in  differentiation  of 
the  spectral  band  with  wavelength  longer  than  twice  the  explo- 
sion depth  (Fuchs,  1966).  Spectral  densities  at  frequencies 
smaller  than  f^  = v^/2H,  where  v^  is  the  average  overburden 
velocity  and  H the  explosion  depth-of-burial , should  decrease 
with  frequency.  For  typical  values  of  v and  H,  f is  near 
1 Hz;  slightly  lower  values  for  very  large  yield  bombs  and 
higher  values  for  very  low-yield  bombs.  While  similar  pP 
cancellation  can  occur  for  earthquakes,  that  portion  of  the 
spectrum  affected  is  generally  outside  the  teleseismically 
recorded  short-period  band  owing  to  the  greater  focal  depths 
of  most  shallow  earthquakes. 

Several  authors  (Wyss,  et  al^,  1971;  Molr.ar,  1971)  have 
pointed  out  that  the  low  frequency  (f  < 1 Hz)  degradation  in 
observed  explosion  spectra  is  usually  faster  than  f,  more  like 
f . This  has  been  attributed  to  the  combined  effects  of  a 
peaked  source  funrtion  and  depth-of-burial. 

Finally,  the  higher  corner/peak  frequencies  for  explo- 
sions as  compared  to  corner  frequencies  for  earthquakes  of 
comparable  is  thought  to  be  a result  of  the  smaller  source 
dimensions  of  explosions  (Wyss  and  Brune,  1968;  Hasegawa, 

1973;  Wyss,  et  al . , 1971).  Rather  dramatic  examples  of  this 
are  given  in  the  paper  by  Wyss,  et  al^  (1971),  where  four 
Aleutian  earthquakes  with  magnitudes  similar  to  the  explosion 
^^row  were  observed  to  have  source  dimensions  at  least  an 
order  of  magnitude  larger  than  those  computed  for  Milrow. 

The  spectral  characteristics  that  were  discussed  in  the 
preceding  paragraphs  for  earthquakes  and  explosions  are  dis- 
played in  Figure  6.1.  This  figure  shows  theoretical  ampli- 
tude spectra  for  an  explosion  and  an  earthquake  as  viewed  at 
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Figure  6.1.  Event  spectra  at  base  of  crust  in  the  source 
region. 
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the  base  of  the  crust  in  the  source  region.  These  two 
events  have  equivalent  body  wave  magnitudes  at  approximately 
1 Hz.  The  effects  of  the  free  surface  reflection  and  crustal 
reverberations  are  included  in  these  spectra. 

The  explosion  spectrum  was  computed  for  a 50  kt  device 
at  a depth-of-burial  less  than  1 kilometer.  The  source 
spectrum  was  comp  i ted  using  a finite  difference  computer 
code  to  model  the  nonlinear  processes  and  elastic  wave  tech- 
niques to  model  the  propagation  of  the  resultant  stress  waves 
to  the  base  of  the  crust  (Bache,  et  al. , 1975  ). 

The  earthquake  spectrum  is  representative  of  a verti- 
cal strike-slip  event  with  a fault  length  of  2 kilometers  and 
a focal  depth  of  15  kilometers.  The  earthquake  source  model 
used  in  this  calculation  (see  Section  IV  of  this  report)  is 
that  of  Archambeau  (1968,  1972)  and  Minster  (1973). 

As  mentioned  previously,  the  VFM  discriminant  is  based 
on  a comparison  of  spectral  magnitude  estimates  made  at  a 
low  frequency  (e.g.,  f ^ 0.4  to  0.5  Hz)  and  a relatively  high 
frequency  (e.g.,  f ^ 2 to  3 Hz).  Examination  of  the  event 
spectra  in  Figure  6.1  indicates  that  the  application  of  the 
VFM  procedure  in  chese  frequency  bands  would  clearly  dis- 
criminate these  two  events.  While  the  holes  or  nulls  in  the 
earthquake  spectrum,  which  are  due  to  reverberations  in  the 
crust,  could  certainly  hamper  discrimination  for  certain 
combinations  of  frequencies,  with  the  use  of  several  fre- 
quencies and  a larger  band  width  between  the  VFM  estimates 
this  potential  problem  can  be  avoided. 

Propagation  of  both  events  through  the  same  upper 
mantle  and  receiver  structure  will  alter  both  spectra  in  the 
same  way;  however,  if  other  stations  or  distances  are  used 
in  combination  for  VFM  discrimination,  then  such  path  effects 
must  be  carefully  taken  into  account. 


6.3 


NARROW- BAND  FILTERING 


Determination  of  the  VFM  discussed  above  is  accomp- 
lished with  the  use  of  a series  cf  narrow-band  filters. 

These  filters  form  the  basic  analytic  portion  of  the  MARS 
signal  processing  code  described  in  cbe  previous  section  of 
this  report.  The  signal  is  filtered  in  the  frequency  domain 
by  multiplication  with  a cusp-type  narrow-band  filter  of  the 
form: 
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The  filter  is  shown  in  Figure  6.2  and  was  chosen  to 
minimize  the  time-domain  ripple  as  discussed  in  Bache,  et  al. 
(197t'  ) . The  width  at  one-half  maximum  amplitude  is  designated 
Af. 


Once  a signal  has  been  narrow-band  filtered,  MARS  then 
corrects  for  the  instrument  response.  If  the  resulting  complex 
spectrum  were  inversely  Fourier  transformed  into  the  time  do- 
main; it  would  represent  a modulated  carrier  wave  of  the 
form: 


x(t)  = A (t)  cos  (uj  t + <t ) , w >0. 

c c 

However,  we  compute  the  Hilbert  transform  of  x(t)  by  multiply- 
ing the  complex  spectrum  by  -i  sgn(w)  and  then  transforming 
into  the  time  domain  by  an  inverse  transform  (see  Savino,  et 
a^-  • > 1975,  for  the  details).  Using  the  signal  and  its 
Hilbert  transform,  the  envelope  and  instantaneous  frequency 
are  produced. 
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Figure  6.2.  Narrow-band  filter  used  in  MARS. 


r 


Corresponding  to  each  filter  center  frequency  and 
envelope  maximum  (Amax)  a seismic  body  wave  magnitude  is 
determined.  The  procedure  followed  in  MARS  is  to  define  the 
oody  wave  magnitude  as 

Vfc)  = lo9  <Amax  * fj  + B (A) 


where 


B(A)  = standard  or  conventional  distance  correction 
factor. 

Each  signal  envelope  maximum  is  corrected  for  noise  by  sub- 
tracting the  maximum  amplitude  of  an  envelope  function,  com- 
puted in  exactly  the  same  manner  as  for  the  signal  envelope, 

based  on  a noise  segment  immediately  preceding  the  signal 
onset. 

Examples  of  narrow-band  filtered  outputs  for  a presumed 
Eurasian  explosion  recorded  by  the  Oyer  subarray  in  Norway 
are  shown  in  Figure  6.3.  In  this  figure  the  unfiltered  time 
series  for  the  explosion  is  shown  in  the  top  left  frame.  The 
remaining  frames  show  the  outputs  of  five  narrow-band  filters 
of  increasing  center  frequency  (f,  = 0.3  to  6.0  Hz)  going 
from  top  to  bottom  and  left  to  right.  The  emergence  of  the 
signal  from  the  noise  at  f,  > l.o  Hz  is  quite  remarkable. 

The  persistence  of  high  frequency  signal  energy  is  likely 
due  to  propagation  along  a high-Q  path. 

As  defined  in  the  previous  paragraph,  the  m,  (f  ) esti- 
mates for  this  event  would  be  based  on  the  maximum  amplitudes 
of  the  filter  envelopes.  As  an  example,  the  m,  (3.0  Hz) 
estimate  would  be  based  on  the  amplitude  of  the  3.0  Hz  fil- 
ter envelope  at  about  27  seconds  into  the  filtered  seismo- 
gram. This  procedure  for  estimating  ^(f..)  will  be  applied 
to  observed  and  synthetic  waveforms  in  the  following  sub- 
sections of  this  report. 
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Figure  6.3.  Examples  of  increasing  signal-to-noise  ratio  for 
a presumed  explosion  (top  left  frame)  recorded 
at  Norway  with  successive  application  of  high 
frequency  narrow-band  filters.  Arrow  at  top 
denotes  approximate  arrival  times  of  the  P— wave 
on  the  unfiltered  trace. 


6 * 4 APPLICATION  of  the  vfm  technique  to  observed  events 

6,4,1  Shallow  Events  Recorded  at  LASA  and  Norway 

The  data  base  employed  in  the  first  test  of  the  VFM 
discriminant  consists  of  LASA  short-period  recordings  of 
P-wave  trains  from  34  presumed  explosions  and  156  earth- 
quakes (Lacoss,  1969).  Thirty  of  the  presumed  explosions 
originated  within  the  mainland  USSR,  two  in  Novaya  Zemlya, 
one  (Longshot)  in  the  Aleutians  and  one  in  the  Sahara  Desert. 
The  earthquakes  are  distributed  along  the  Alpine  seismic 
belt,  the  Kuril-Kamchatka  arc  and  the  Arctic  Ocean.  The  epi- 
central  distances  of  these  events  range  from  45°  to  about 

100  with  more  than  half  of  the  events  between  65°  and  85° 
from  LASA. 

_ Figure  6.4  is  a plot  o::  spectral  magnitude  estimates, 

rob(fc),  at  a low  frequency  (f,  = 0.45  Hz)  versus  a high 
frequency  (fc  - 2.25  Hz)  for  all  34  presumed  explosions  and 
those  earthquakes  in  the  data  base  that  were  reported  in  the 
Monthly  Listing  f Events,  published  by  the  USGS  (formerly 
NOAA) , as  having  snallow  focal  depths,  h < 70  km.  This 
figure  clearly  demonstrates  the  enriched  high  frequency 
content  of  the  explosion  body  wave  spectra  as  compared  to 
the  earthquake  spectra.  For  instance,  for  a given  value  of 
”b  (0,45  Hz)  the  explosions  exhibit  ir^  (2.25  Hz)  values  that 
are  typically  0.6  to  1.0  unit  larger  than  the  m.  (2.25  Hz) 
values  for  earthquakes.  ® 

The  high  degree  of  discrimination  of  earthquakes  from 
explosions  evident  in  Figure  6.4  is  especially  significant 
m view  of  the  non-regionalization  of  the  event  population. 

The  variety  of  tectonic  settings  of  this  event  population 
ranges  from  relatively  stable  shield  regions  to  seismically 
active  oceanic  arc  systems.  An  indication  that  discrimination 
could  be  further  enhanced  by  regionalizing  the  event  popula- 
tion comes  from  the  fact  that  for  ir^  (2.25  Hz)  >4.0  the  two 
presumed  explosions  (#35  and  #138)  plotting  closest  to  the 
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Spectral  magnitudes,  mjj,  computed  at  0.45  Hz 
and  2.25  Hz.  The  presumed  explosions  numbered 
35  and  138  occurred  at  Novaya  Zemlya. 


earthquake  population  occurred  at  Novaya  Zemlya . In  addi- 
tion, several  of  the  earthquakes  that  plot  closest  to  the 
explosion  population  occurred  along  the  Kurile-Kamchatka 
arc  or  at  locations  far  removed  from  the  explosion  epicenters. 

A subset  of  the  shallow  event  ensemble  recorded  by 
LASA  was  also  recorded  by  the  limited  Oyer  array  at  Norway. 

The  VFM  technique  was  applied  to  this  small  subset  of  events 
and  the  results  are  shown  in  Figure  6.5.  Analogous  to  the 
results  for  LASA,  there  is  discrimination  of  events  over 
most  of  the  magnitude  range  except  for  the  smallest  magni- 
tude explosion. 

The  apparent  bending  of  the  explosion  population  into 
the  earthquake  population  at  m^  (2.25)  < 3.5  at  LASA  and 
m^  (5.0)  < 4.0  at  Norway  is  mainly  the  result  of  microseis- 
mic  noise  enhancing  the  low  frequency  ir^  estimates  for  the 
small  (signal-to-noise  sense)  explosions.  Correction  for 
this  noise  at  Norway,  results  in  complete  separation  of 
earthquakes  and  explosion  populations  (Figure  6.6).  There 
was  insufficient  noise  data  available  for  LASA  to  make  simi- 
lar corrections  for  noise  to  that  data  set. 

For  a given  recording  site  the  limitation  of  the  low 
frequency  end  of  the  seismic  bandpass  of  the  VFM  discriminant 
is  controlled  by  the  microseismic  noise  field  at  that  station. 
At  the  high  frequency  end  of  the  seismic  bandpass  there  are 
three  factors  that  limit  the  VFM  discriminant  to  small  magni- 
tude teleseismic  events.  These  three  factors  are  source 
spectrum,  anelastic  attenuation  and  high  frequency  background 
noise. 


Extensive  study  of  various  VFM  discriminant  parameters 
indicates  that  optimum  discrimination  would  result  for  spec- 
tral magnitudes  computed  at  frequencies  (fc)  separated  by  a 
decade  or  more  such  as  observed  at  Norway.  However,  noise 
and  anelastic  atcenuation  of  the  earth  combine  to  constrain 
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Spectral  magnitude  estimates  at  fc  = 0.6  Hz 

andjfS  = 5.0  Hz  for  an  event  population  re- 
corded at  the  Oyer  array  in  Norway.  The 
arrows  attached  to  several  of  the  closed 
ci^c^es  indicate  that  the  low  frequency  (O.i 
Hz)  mb  estimates  are  at  the  prevailing  noise 
level  as  determined  from  recording  segments 
immediately  preceding  each  signal  onset. 
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Figure  6.6.  Maximum  filter  amplitudes  with  noise  correction  for  the  same  event 
population  plotted  in  Figure  6.5. 


the  usable  bandwidth  of  frequencies  in  the  case  of  LASA  data 
set  to  approximately  one-half  a decade. 

6.4.  2 VFM  for  Deep  Events  (h  > 70  km) 


At  LASA,  twenty  events  in  the  depth  range  of  80  km 
to  580  km  were  included  in  the  data  base  and  at  Norway 
eleven  events  with  depths  greater  than  70  km  were  included. 
Figures  6.7  and  6.8  shew  the  results  of  the  application  of 
VFM  to  these  data  as  compared  to  the  results  obtained  ior 
shallow  events  (contoured  solid  and  dashed  lines). 

Deep  earthquakes  often  fail  to  separate  from  explosion 
populations  when  examined  with  any  of  the  discriminants, 
short-  and/or  long-period,  proposed  to  date.  A number  of 
reasons  have  been  proposed  for  the  failure  of  the  discriminants 
when  applied  to  many  deep  earthquakes.  These  include  the 
hypotheses  that  the  deep  earthquakes  are  characterized  by 
impulsive  source  time  functions  and/or  small  source  dimen- 
sions or  that  the  source-receiver  paths  are  characterized  by 
high  Q. 

While  some  of  the  deep  earthquakes  overlap  iro  the 
explosion  populations,  they  are  really  not  as  trouh  -some  as 
they  might  seem  since  these  events  can  be  identified  as 
naturally  occurring  earthquakes  on  the  basis  of  other  infor- 
mation (e.g.,  hypocentral  location,  using  P,  pP,  and/or 
sP  travel  times) . 

A rather  dramatic  example  of  the  short-period  spectral 
similarities  between  an  explosion  and  a deep  earthquake  is 
shown  in  Figure  6.9.  Approximately  80  seconds  of  an  un- 
f liter ed  explosion  time  series  and  the  output  from  three 
narrow-band  filters  (fc  = 0.3,  1.0  and  6.0  Hz)  applied  to 
the  time  series  are  shown  on  the  left-hand  side  of  the 
figure.  A similar  sequence  of  pictures  is  shown  on  the  right- 
hand  side  for  a deep  earthquake  (h  = 511  km) . Note  the 
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Figure  6.7.  Spectral  magnitude  estimates  of  deep  earthquakes 
recorded  at  LASA.  The  shallow  earthquake  and 
explosion  populations  plotted  in  Figure  6.4 
are  contoured  in  this  figure. 
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Figure  6.9.  Narrow-band  filter  outputs  at  three  different 
center  frequencies  (fc  = 0.3,  1.0  and  6.0  Hz) 
for  an  explosion  (left-hand  side)  and  a deep 
earthquake  (right-hand  side) . 
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enhancement  of  signal-to-noise  ratio  for  the  explosion  and 
earthquake  signals  at  the  higher  frequencies  (f  > 1.0  Hz). 

6*4*3  North  American  Events  Recorded  in  Canada 

Data  in  digital  format  for  a large  population  of  North 
American  earthquakes  and  underground  explosions  recorded  at 
the  19-element  short-period  Yellowknife  array  in  Canada 
(Manchee,  1972)  were  recently  acquired.  These  data  consist 
of  vertical-component  recordings  of  P-wave  trains  (preceded 
by  at  least  25  seconds  of  background  noise)  from  approximately 
40  earthquakes  occurring  in  the  Gulf  of  California  and  49 
underground  explosions  detonated  at  the  Nevada  Test  Site 
(NTS) . While  signals  recorded  on  all  19  array  elements  were 
supplied  for  most  of  the  events  in  the  data  base,  our  first 
test  for  discrimination  with  the  VFM  technique  was  performed 
on  single-element  (sensor  B5  located  at  the  crossover  point 
of  the  L-shaped  array  arms)  recordings  only. 

Each  of  the  time  series  tested  for  discrimination  in 
this  study  was  first  divided  into  a noise  sample  and  a 
signal  sample.  In  the  case  of  the  explosions,  both  noise 
and  signal  windows  were  chosen  to  be  24  seconds  long.  For 
the  earthquakes,  a 24  second  noise  window  and  signal  windows 
up  to  66  seconds  in  duration  were  examined.  The  longer 
duration  of  the  earthquake  signal  windows  was  dictated  be- 
cause of  the  difficulty  in  identifying  the  exact  signal 
onset  time. 

Each  of  the  signal  time  series  was  filtered  by  twenty 
narrow  band  filters  with  center  frequencies  ranging  from 
0.25  Hz  to  4.0  Hz.  The  frequency  range  was  selected  after 
examining  spectra  of  ten  different  earthquakes  and  explo- 
sions. Magnitude  estimates,  n^ff),  based  on  a narrow  band 
filter  with  a relatively  low  center  frequency,  f , were 
then  plotted  versus  estimates  based  on  high  frequency 
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filters.  This  was  repeated  for  various  combinations  of 
center  frequencies  and  all  of  the  plots  were  then  examined 
for  discrimination. 

Figure  6.10  shows  results  for  the  North  American  data 
base.  The  filter  center  frequencies  that  have  yielded  the 
best  discrimination  results  are,  as  indicated,  0.425  Hz  and 
2.5  Hz.  In  this  figure  the  Gulf  of  California  earthquakes 
are  denoted  by  the  open  circles,  the  NTS  explosions  by  the 
X s.  In  general  most  of  the  earthquakes  on  this  figure, 
which  are  reported  in  the  Preliminary  Determination  of  Epi- 
centers to  be  shallow  focus,  separate  from  the  explosion 
population.  While  the  results  in  Figure  6.10  are  not  nearly 
as  impressive  in  terms  of  complete  separation  of  eve nh 
populations  (e.g.,  as  in  Figure  6.4  for  Eurasian  events  at 
LASA) , a very  important  and  positive  point  about  the  dis- 
crimination test  in  Figure  6.10  is  the  complete  separation 
of  the  very  small  magnitude  events. 


6 • 5 MULTIPLE  EXPLOSION  SCENARIO 

The  VFM  technique  is  especially  suited  for  identifica- 
tion and  discrimination  of  multiple  explosion  sequences  that 
are  designed  to  appear  earthquake-like  in  terms  of  the  con- 
ventional (Mg-n^,  depth  of  focus,  complexity,  first  motion) 
discriminants.  A multiple  event  scenario,  somewhat  similar 
to  one  proposed  by  Kolar  and  Pruvost  (1S75) , was  devised  by 
superposing  eight  scaled  seismograms  of  a presumed  Kazakh 
explosion  recorded  at  LASA.  This  explosion  signature  is  shown 
on  the  bottom  center  of  Figure  6.11. 

The  array  of  explosions  and  their  relative  yields  were 
designed  to  produce  earthquake-like  seismograms  over  a wide 
range  of  azimuths.  The  particular  array  configuration  (spac- 
ing and  firing  order)  is  indicated  in  the  center  of  the  figure. 
Each  seismogram  comprising  the  multiple  event  was  delayed  in 
time  relative  to  the  first  end  scaled  in  amplitude.  The 
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Figure  6.11.  Primary  explosion  signature  (bottom-center)  used 
to  make  composite  seismograms  at  five  different 
azimuths  from  the  array  of  eight  explosions. 

The  firing  order  and  spacing  of  the  explosions 
are  indicated  in  the  center. 
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explosion  delay  times  between  sources  are:  0.0,  0.35,  1.2 

2.9,  4.7,  5.7,  6.0  and  6.3  seconds,  respectively.  The  <m-' 
plltude  scaling  for  all  eight  explosions  is  1,  3.1,  5,  10 
10,  20,  15  and  12.5.  The  largest  explosion  in  the  set  is' 

e sixth  event  and  is  scaled  to  give  the  same  ground  motion 
as  tile  primary  signal. 

The  composite  seismograms  are  shown  in  Figure  6.11 
for  five  different  azimuths  (1-5)  with  respect  to  the  shot 
array.  The  first  point  to  be  noted  about  these  composite 
seismograms  is  that  with  the  addition  of  noise  at  the  begin- 
ning of  each  seismogram,  the  first  motion  at  azimuths  1,  3 
and  5 would  most  likely  be  picked  as  rarefactions.  Secondly, 
the  complexity  of  each  composite  signal  has  been  significant!, 
increased  over  that  of  the  primary  signal. 

Two  factors  work  together  to  make  the  multiple  event 
scenario  appear  earthquake-like  on  the  M -m.  basis.  First 
the  analyst  making  amplitude  measurements  in  the  conven- 
tional manner  for  the  determination  of  - would  undoubtedly 
Pick  amplitudes  corresponding  to  the  earlier  smaller  explo- 
sion in  the  sequence.  On  the  other  hand,  the  M measure- 
ments would  be  based  on  the  superposed  surface  waves  from  the 
ree  largest  explosions  occurring  later  in  the  sequence. 

The  results  would  be  to  reduce  the  at  and  enhance  the  M 
estimates  and  place  the  composite  event  in  the  earthquake 
population  on  an  Mg-m^  plot. 

However,  application  of  the  VFM  technique  results  in 
complete  discrimination  of  the  multiple  explosion  scenario. 

" 19Ure  6-12  the  VFM  estimates  determined  for  the  five 
azimuths  are  indicated  by  the  X's.  A most  significant  re- 
sult is  that  the  VFM  estimates  of  the  multiple  explosion 
sequence  at  the  five  azimuths  cluster  around  the  estimates 

° e primary  S1«nal  used  to  form  these  composite  signals 
(the  closed  circle  immediately  to  the  right  of  x1.  This 

means  that  the  VFM  technique  has  based  the  magnitude 
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Figure  6.12.  Spectral  magnitudes  for  same  event  population 
and  filter  parameters  as  in  Figure  6.4  with 
estimates  for  multiple  explosions  (X's)  showing 
complete  discrimination. 
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estimates  on  the  largest  amplitude  arrival  in  the  wavetrain, 
corresponding  to  the  largest  yield  explosion  in  the  sequence. 
This  is  an  important  result  for  yield  determination  of  both 
single  and  multiple  explosions. 


6,6  application  of  vfm  to  theoretical  seismograms 

In  this  section  we  compute  VFM  values  for  a suite  of 
synthetic  earthquake  and  explosion  seismograms.  Our  purpose 
in  doing  this  is  twofold.  First,  theoretical  seismograms 
provide  a controlled  data  set  for  testing  and  improving  a 
data  analysis  technique,  such  as  the  VFM.  With  theoretical 
seismograms  the  controlling  parameters  and  their  effect  are 
thoroughly  understood.  Second,  we  want  to  compare  the  theo- 
retical VFM  values  to  the  observed  values  to  see  how  well 
our  theoretical  seismograms  match  actual  data  with  regard  to 
this  spectral  measure. 

This  section  is  divided  into  four  parts.  We  first 
briefly  describe  the  computation  of  the  synthetic  seismograms 
used  in  the  study.  Then  the  VFM  estimates  are  computed,  first 
for  explosions,  then  for  earthquakes.  Finally,  the  theore- 
tical VFM  values  are  compared  to  those  from  observations. 

6«6.1  Calculation  of  Synthetic  Seismograms 

The  techniques  used  to  compute  short  period  body  wave 
theoretical  seismograms  have  been  developed  under  this  con— 
^rac^  anc*  have  been  described  in  previous  semi-annual  reports 
(Bache,  et  ah,  1974,  1975;  Savino,  et  al.,  1975).  The  tech- 
nique is  outlined  in  Chapter  IV  of  this  report  and  a number 
of  examples  are  given  there.  Here  we  shall  rely  upon  the 
discussion  in  Chapter  IV  and  simply  point  out  features 
peculiar  to  the  synthetic  seismograms  for  the  VFM  study. 

Synthetic  seismograms  were  computed  for  five  dif- 
ferent events  which  are  intended  to  be  typical  of  NTS.  The 
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characteristics  of  these  five  events  are  summarized  in 
Table  6.1. 

Table  6.1.  Characteristics  of  Explosion  Events  for 
Synthetic  Seismograms 


Event 

Test  Area 

Source 

(#) 

Yield 

(kt) 

Burial  Depth 
(km) 

1 

Yucca  Flat 

Dry  Tuff 

(134) 

10 

0.35 

2 

Pahute  Mesa 

Dry  Tuff 

(134) 

50 

0.42 

3 

Pahute  Mesa 

Wet  Tuff 

(131) 

400 

0.875 

4 

Pahute  Mesa 

Rhyolite 

(125) 

1000 

1.20 

5 

Pahute  Mesa 

Rhyolite 

(125) 

1500 

1.35 

The  coupling  of  explosion  energy  into  elastic  waves 
varies  with  the  characteristics  of  the  material  type.  In  our 
computations  the  equivalent  elastic  source  is  specified  by 
a reduced  displacement  potential  which  is  computed  using 
finite  difference  methods.  The  procedure  for  computing  the 
source  has  been  described  by  Cherry,  et  al^(l974).  The 
sources  134,  131  and  125  were  computed  for  NTS  dry  tuff,  wet 
tuff  and  rhyolite,  respectively,  and  are  discussed  by  Bache, 
et  al^  (1975) . 

For  the  synthetic  seismogram  calculations  two  source 
region  crustal  structures  were  used,  one  for  Pahute  Mesa  and 
one  for  Yucca  Flat.  These  are  tabulated  in  Tables  6.2  and 
6.3. 

Except  for  the  source  specification  in  terms  of  a 
reduced  displacement  potential  and  the  velocity  profile  in 
the  source  region,  the  synthetic  seismogram  calculations 
are  identical  to  those  for  earthquakes  described  in  Chapter 
IV.  Once  again,  the  upper  mantle  model  is  HWNE  (Helmberger 
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Table  6.2.  Yucca  Plat  -rust  for  Theoretical  Seismograms 

P Wave  S Wave 

Thickness  Velocity  Velocity  Density 
.yanj — (fan) (fan/sec)  (km/<?»rO  (gm/cm*) 


Layer  (km) 
1 0.3 


(fan) 

0.3 


(km/ sec)  (km/ sec) 
1.585  0.9 


1.65 


2 

0.375 

0.075 

2.286 

1.212 

1.73 

7 

0.6 

0.225 

2.13 

0.99 

1.7 

4 

0.955 

0.325 

2.5 

1.43 

1.9 

5 

1.175 

0.25 

4.4 

2.75 

2.2 

6 

2.3 

1.125 

5.0 

3.1 

2.4 

7 

4.3 

2.0 

5.5 

3.2 

2.6 

8 

20.0  15.7 

6.0 

3.5 

2.8 

Table 

6.3. 

Pahute  Mesa  Crust  for  Theoretical  Seismograms 

Layer 

Depth 

P Wave 

Thickness  Velocity 

S Wave 

Velocity  Density  Material 

1 

0.256 

0.256 

2.74 

1.3 

2.0 

2 

0.36 

0.104 

3.6 

1.9 

2.3 

3 

0.65 

0.29 

2.13 

1.1 

2.1 

Dry  Tuff 

4 

0.8 

0.15 

2.5 

1.56 

1.85 

Wet  Tuff 

5 

6 
7 

1.4 

2.1 

6.0 

0.6 

0.7 

3.9 

4.2 

4.3 
4.7 

2.62 

2.4 

2.6 

2.4 

2.6 

2.6 

Wet  Rhyolite 
Wet  Rhyolite 

8 

12.0 

6.0 

5.4 

2.7 

2.7 

9 

20.0 

8.0 

6.0 

3.5 

2.8 
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and  Wiggins,  1971)  and  an  average  continental  crust  (Table 
r 4.2)  is  used  at  the  receiver.  The  T/Q  is  1.05,  a value 

found  to  be  appropriate  for  NTS  northeast  profiles. 

At  a single  epicentral  distance  seismograms  for 
events  1 - 5 of  Table  6.1  will  illustrate  the  effect  of 
variations  in  the  source  function  (reduced  displacement 
potential),  the  yield  and  the  burial  depth  or  P-pP  lag 
time.  We  also  would  like  to  see  the  effect  of  differences 
in  the  upper  mantle  travel  path  so  seismograms  were  com- 
puted for  two  distances;  3735  and  2360  km.  The  synthetic 
seismograms  are  shown  in  Figure  6.13.  At  3735  km  the  upper 
mantle  response  is  beyond  the  triplication  range  and  the 
seismograms  are  quite  simple.  The  response  at  2360  km  is 
within  the  triplication  range  and  the  seismograms  are  con- 
siderably more  complex. 

The  computation  of  synthetic  seismograms  for  earth- 
quakes using  the  Archambeau/Minster  source  model  was  de- 
scribed in  detail  in  Chapter  IV.  In  that  ‘chapter  the 
scaling  of  seismogram  amplitudes  with  the  source  parameters 
is  discussed.  Using  these  techniques  and  precisely  the  same 
earth  models  used  for  the  seismograms  of  Chapter  IV,  a suite 
of  theoretical  seismograms  were  computed  in  which  the  para- 
meters varied  include  the  depth  of  focus,  the  fault  length, 
the  fault  orientation  and  the  epicentral  distance.  We  also 
should  mention  that  the  form  of  the  Archambeau/Minster 
source  used  for  these  calculations  is  not  identical  to  that 
of  Chapter  IV,  but  that  the  spherical  volume  of  vanishing 
rigidity  both  expands  and  propagates  (in  the  model  used  in 
Chapter  IV  the  spherical  volume  expands  from  a fixed  point 
on  the  diameter) . For  our  purposes  the  differences  between 
the  models  are  quite  small  and  the  scaling  behavior  described 
in  Chapter  IV  is  valid  for  both.  The  material  properties, 

rupture  velocity,  etc.,  are  the  same  as  for  the  earthquakes 
of  Chapter  IV. 
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As  in  Chapter  IV  we  compute  synthetic  seismograms 
for  the  five  focal  depths  10,  15,  25,  35  and  45  km.  Earth- 
quake sources  were  computed  at  three  different  fault  lengths 
In  the  material  at  35  - 45  km  depth  these  lengths  were  10, 

5 and  2.5  km.  Then  according  to  the  scaling  properties  of 
the  model  discussed  in  Section  4.2,  the  fault  lengths  are 
8.4,  4.2,  2.1  kn  at  25  km  depth  and  7.6,  3.8  and  1.9  km  at 
the  10  and  15  km  depths.  All  seismograms  were  scaled  to  a 
100  bar  stress  drop. 

In  Figures  6.14  and  6.15  we  show  typical  examples  of 
the  synthetic  seismograms  for  this  VFM  study.  The  only  dif- 
ference between  the  seismograms  of  Figures  6.14  and  6.15  is 
in  the  upper  mantle  response  which  in  one  case  is  for  a 
range  of  2500  km,  within  the  triplications,  and  for  the 
other  case  is  for  a range  of  4000  km.  The  source  orienta- 
tion is  right-lateral  strike-slip  and  the  azimuth  is  30° 
counterclockwise  from  the  strike  (see  Figure  4.2).  As  for 
the  explosions,  the  seismograms  at  the  nearer  range  are  con- 
siderably more  complicated. 

Seismograms  for  the  three  basic  sources  described 
above  were  also  computed  for  two  other  source  orientations 
in  addition  to  the  strike-slip.  These  orientations  are 
dip-slip  and  45°  thrust.  For  both,  the  rupture  initiates 
at  the  origin  and  breaks  downward.  For  the  dip-slip  events 
the  fault  plane  shown  in  Figure  4.2  is  rotated  90°  about 
the  Y_-axis.  That  is,  the  material  on  the  +Y_  side  of  the 
XG  - ZG  plane  is  downthrown.  For  the  thrust  events  the 
dip-slip  fault  plane  is  rotated  45°  clockwise  about  Xq. 

Then  the  hanging  wall  slides  upward  over  the  footwall. 

6.6.2  VFM  Values  for  Synthetic  Explosion  Seismograms 


Fox  the  synthetic  seismograms  shown  in  Figure  6.13, 
VFM  values  were  computed  at  a number  of  frequencies.  For 


these  events  we  plot  (0.45  Hz)  versus  m,  (2.25  Hz)  in 
Figure  6.16.  “ 


Recall  that  the  events  for  which  the  synthetic  seis- 
mograms were  computed  varied  in  emplacement  material,  depth 
of  burial  and  yield.  Also,  seismograms  were  computed  at 
two  epicentral  distances.  The  effect  of  the  distance  varia- 
tion is  clear  and  consistent  with  the  (f ) about  the  same 
for  every  event.  The  lower  frequency  ra^f)  values  are  much 

less  affected  by  the  introduction  of  the  upper  mantle  trip- 
lications. 


The  variation  of  the  mb(f)  values  with  the  source 
parameters  is  strongly  dependent  on  frequency.  As  one  might 
expect,  the  low  frequency  values  (it^  (0.45  Hz))  are,  pro- 
portional to  the  yield.*  The  high  frequency  values  on  the 
other  hand  seem  to  be  almost  independent  of  yield.  This  is 
because  the  high  frequency  (2.25  Hz)  is  beyond  the  corner 
frequency  for  these  events. 


6 * 6 * 3 Values  for  Synthetic  Earthquake  Seismograms 

For  the  earthquake  synthetic  seismograms  the  parameters 
varied  include  the  focal  depth,  fault  length,  fault  orienta- 
tion and  epicentral  distance.  The  stress  drop  is  held  constant 
and  the  other  fault  parameters  obey  the  similarity  rules  set 
forth  in  Chapter  IV. 


First,  let  us  examine  the  effect  of  fault  orientation 
on  the  VFM  estimates.  The  basic  event  is  that  which  has  a 

T~ 

Actually  the  proportionality  is  to  the  product  of  the  local 
P wave  veiocity  and  the  amplitude  of  the  reduced  velocity 
the_£ominant  frequency  as  is  shown  by  Bached 

,^75  The  source  function  amplitude  is,  of  course, 
strongly  dependent  on  the  explosion  yield. 
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Figure  6.16.  Theoretical  NTS  explosion  VFM  estimates  for  two 
distances  (21°,  36°). 
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fault  length  of  5.0  km  at  the  45  km  depth.  The  n^(f)  were 
computed  from  synthetic  seismograms  at  a single  azimuth  (30°) 
and  distance  (4000  km).  For  the  strike-slip  orientation  the 
seismograms  were  shown  in  Figure  6.14. 

In  Figure  6.17  we  plot  (0.45)  versus  5^  (2.25)  for 
the  fifteen  synthetic  seismograms.  In  terms  of  our  discus- 
sion of  the  scaling  of  teleseismic  body  wave  amplitudes  in 
Section  4.6,  all  these  events  have  the  same  "source  strength". 
Then  the  simple  scaling  relationship  (4.13)  indicates  that, 
fixing  the  radiation  pattern  factor,  the  amplitude  should 
scale  with  a or  L since  L/a  is  fixed.  Except  for  the 
strike-slip  orientation  we  see  that  the  VFM  values  are 
poorly  explained  by  this  scaling.  This  is  apparently  be- 
cause the  radiation  pattern  factor  is  strongly  varying  from 
record-to-record,  even  for  a fixed  fault  orientation. 

In  Figure  6.1*  we  plot  the  VFM  values  for  the  three 
different  source  strength  events  under  study.  Once  again 
we  show  the  values  for  all  five  focal  depths.  The  source 
orientation  is  strike-slip  and  the  azimuth  and  distance  are 
again  30°  and  4000  km.  Within  each  group  the  variation 
with  depth  (and  length)  is  about  the  same  but  the  major 
change  is  due  to  the  change  in  source  strength. 

Finally,  in  Figure  6.19  we  select  one  of  the  sources, 
the  strike-slip  earthquake  with  L = 5.0,  4.2,  3.8  km,  and 
show  the  effect  of  changing  epicentral  distance.  That  is, 
the  seismograms  are  those  of  Figures  6.14  and  6.15.  As  with 
the  explosions,  we  see  that  the  effect  of  the  upper  mantle 
response  variation  is  quite  consistent  from  event  to  event. 

In  summary , we  have  seen  that  the  variation  of  VFM 
with  the  source  parameters  is  quite  complex  and  not  easily 
explained  with  simple  scaling  relationships.  The  data  for 
superficially  similar  events  scatters  over  a range  of  more 
than  0.5  VFM  units  at  both  high  and  low  frequencies.  In  the 
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Figure  6.17.  Plot  of  VFM  for  five  events  at  three  different 
fault  orientations.  The  distance  and  azimuth 
are  4000  km  and  30°  for  all  seismograms. 
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Figure  6.18.  Theoretical  VFM  estimates  for  three  strike- 
slip  earthquakes. 
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Figure  6.19.  Theoretical  VFM  estimates  for  the  seismograms 
of  Figures  6.14  and  6.15. 
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next  section  we  consider  the  date  in  a global  sense  and  see 

how  well  the  trends  of  our  theoretical  estimates  agree  with 
observations. 


6-7  COMPARISON  OF  SYNTHETIC  TO  OBSERVED  VFM  nama 

This  section  presents  a comparison  of  theoretical  ex- 
plosion and  earthquake  VFM  estimates  to  selected  observed 
VFM  data  discussed  previously  (see  Section  6.4).  The  ob- 
served data  from  explosions  includes  NTS  explosions  recorded 
at  YKA  (Figure  6.10).  Observed  earthquake  VFM  data  are  from 
the  Gulf  of  California  and  recorded  at  YKA  (Figure  6.10). 


6*7#1  Comparison  of  Explosion  VFM  Data 

In  Figure  6.20  we  show  the  theoretical  explosion  data 
discussed  in  the  previous  section  plotted  with  the  observed 
NTS  explosion  data  recorded  at  YKA.  The  distance  range  for 
the  observed  YKA  data  is  25.2"  to  25.4"  versus  the  distance 
range  of  21"  to  36«  for  the  theoretical  data.  The  agreement 
between  the  observed  and  synthetic  explosion  populations  in 
terms  of  absolute  magnitude  estimates  and  trend  is  quite 
good.  For  instance,  three  of  the  largest  observed  explosions 
were  in  the  150  to  200  kt  yield  range  while  ten  had  reported 
yields  of  10  kt  or  less.  The  majority  of  the  remaining  event 
were  in  the  50  to  100  kt  yield  range.  Further,  note  that  the 
noise  corrections  applied  to  the  observed  data  set  were 
particularly  significant  for  the  lower  yield  bombs.  In 
general  these  corrections  resulted  in  a reduction  of  the  VFM 

estimates  which  accounts,  in  part,  for  the  very  low  magnitude 
estimates  in  Fig-are  6.20. 


In  the  previous  section  of  this  report  we  showed  that 
there  is  a strong  dependence  of  the  explosion  VFM  estimates 
on  depth-of-burial  and  near  source  material  properties.  In 
view  of  the  fact  that  the  explosion  population  in  Figure 
6.20  includes  events  from  several  different  test  areas  at 
NTS  (i.e.,  Pahute  Mesa,  Yucca  Valley,  Rainier  Mesa),  the 
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• Theoretical  Explosions,  A « 21c 
□ Observed  NTS  Explosions  at  YKA,  a 


(2.25  Hz) 


Figure  6.20. 


<;f^aris?n  of  theoretical  to  observed  explosion 
(NTS-YKA)  VFM  estimates. 


amount  of  scatter  of  the  observed  data  is  not  really  surpris- 
ing. The  scatter  could  undoubtedly  be  reduced  by  segregating 
the  explosions  according  to  the  specific  test  areas  and  com- 
paring with  appropriate  synthetic  events.  However,  our  in- 
tent here  was  to  approximate  the  gross  features  of  the  ob- 
served data;  an  objective  which  is  obviously  fulfilled  in 
Figure  6.20. 


6,7,2  Comparison  of  Earthquake  VFM  Data 

The  theoretical  VFM  estimates  for  the  strike-slip 
earthquakes  (0%  90%  0%  at  a distance  of  36“  shown  in 
Figure  6.18  for  three  fault  sizes  (10  * 8.4  km,  5 * 4.2  km 
and  2.5  x 2.1  km)  and  five  depths  (10,  15,  25,  35  and  45  km) 
are  plotted  in  Figure  6.21  along  with  the  Gulf  of  California 
VFM  data.  The  average  distance  of  YKA  from  the  Gulf  of 
California  is  approximately  37 • versus  36*  for  the  theoreti- 
cal data.  The  theoretical  values  are  positioned  within  the 
left  portion  of  the  Gulf  of  California  values. 

The  scatter  of  the  observed  data  is  large.  Thatcher 
and  Brune  (1971)  indicate  that  characteristically  these  earth- 
quakes are  very  shallow  and  are  predominantly  normal  faults. 
Had  we  used  dip-slip  earthquakes  (0%  90%  90%  theoretical 
VFM  values  instead  of  the  strike-slip  type,  the  theoretical 
population  would  be  positioned  near  the  center  of  the  ob- 
served data.  The  theoretical  results  shown  previously  in 
Figure  6.17  confirm  this  statement.  However,  much  more  must 
be  learned  about  the  fault  mechanisms,  stress-drops  and  depths 
of  these  Gulf  of  California  earthquakes  to  explain  the  large 
variance  observed  in  the  VFM  estimates. 

6,7,3  Comparison  of  Explosion-Earthquake  VFM  Data 

Figure  6.22  presents  an  overview  of  the  important  data 
'ncluded  in  this  report.  In  this  figure  the  theoretical  VFM 
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Figure  6.22. 


Comparison  of  theoretical  and  observed 
sion  and  earthquake  VFM  estimates. 


data  for  explosions  and  strike-slip  earthquakes  are  plotted 
corresponding  to  an  epicentral  distance  of  4000  kms.  In 
addition  we  outline  the  Gulf  of  California  population 
(broken  line)  and  the  NTS  explosion  (solid  line)  as  observed 
at  YKA.  We  previously  indicated  that  had  we  used  theoreti- 
cal dip-slip  earthquakes  the  VFM  estimates  would  lie  near 
the  center  of  the  observed  Gulf  of  California  VFM  data. 
Further,  the  observed  NTS-YKA  data  occurs  at  about  25.3° 
distance  and  plotting  the  theoretical  values  at  a distance 
of  21°  on  the  figure  (solid  triangles)  indicates  that  they 
too  lie  near  the  center  of  the  observed  explosion  popula- 
tion. Hence,  the  mean  separation  between  earthquakes  and 
explosions  populations  of  both  theoretical  and  observed 
groups  is  approximately  the  same  (=  0.8  units  of  magnitude 
normal  to  the  population  trend  lines) . 

6.8  SUMMARY  AND  CONCLUSIONS 

Application  of  the  variable  frequency  magnitude  (VFM) 
discriminant  to  North  American  and  Eurasian  event  popula- 
tions yields  complete  separation  between  earthquakes  and  ex- 
plosions down  to  conventional  magnitudes,  m^,  in  the  range 
of  4.0  to  4.5.  Enhancement  of  this  discriminant  can  be 
achieved  by  regionalization  of  the  event  populations  and 
the  application  of  noise  corrections  to  the  VFM  values. 
Obviously,  the  latter  process  becomes  most  important  for 
small  events. 

Application  of  the  VFM  discriminant  to  a multiple 
explosion  scenario  yielded  VFM  values  similar  to  those  of 
single  explosions  for  five  azimuths  sampled  around  the 
linear  array  of  explosions.  All  VFM  values  clustered  arour.d 
that  of  the  primary  signal  employed  and  correspond  to  the 
largest  yield  explosion  in  the  sequence. 

Estimation  of  the  VFM  values  for  theoretical  explo- 
sions and  earthquakes  originating  in  the  same  source  region 
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and  observed  at  similar  ep„central  distances  indicate  the 
following: 

• Explosion  VFM  exhibit  a strong  dependence  on 

depth-of-burial  and  near  source  material  pro- 
perties. 

e For  equal  sized  earthquakes  of  various  orienta- 
tions the  most  significant  variation  in  the  VFM 
estimates  is  associated  with  focal  depth.  Ob- 
servations of  similar  fault  orientations  of 
varying  size  shows  that  depth  changes  result  in 
nonuniform  variations  in  the  VFM  estimates. 

Comparison  of  VFM  values  for  theoretical  and  observed  explo 
srons  and  earthquakes  where  Q's  and  distance  ranges  are 
reasonably  comparable  show  the  following: 

• VFM  estimates  for  NTS  explosions  observed  at 
YKA  compare  well  to  the  .theoretical  VFM 
values  in  terms  of  both  absolute  values  and 
the  trends  of  the  two  populations. 

• Gulf  of  California  VFM  observed  at  YKA  and 
compared  to  strike-slip  earthquakes  did  not; 
compare  well.  However,  theroetical  dip-slip 
earthquake  VFM  values  are  positioned  near 
the  center  of  the  observed  data;  in  agreement 
with  the  observation  that  most  of  the  Gulf  of 
California  earthquakes  were  predominantly  dip- 
slip  events. 

Comparing  all  explosion-earthquake  data  (theoretical 
and  observed)  we  find  that  the  mean  separation  between  earth 
quakes  and  explosions  to  be  approximately  0.8  units  of 
magnitude  units  normal  to  the  population  trends  and  for 
frequencies  of  0.45  and  2.25  Hz. 
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Future  work  should  concentrate  on  a systematic  study 
of  the  Gulf  of  California  earthquake  parameters  (i.e 
depth,  stress  drops  and  fault  orientations)  to  determine  the 
causes  of  the  large  variance  seen  in  the  VFM  estimates. 
Further,  we  have  investigated  only  a few  theoretical  explo- 
sion and  earthquake  parameter  changes  and  reported  their 
effects  on  the  VFM  estimate  in  a qualitative  way.  We  need 
to  know  what  parameters  are  particularly  sensitive  in  the 
estimation  of  VFM  to  determine  just  how  much  confidence  can 
be  placed  on  this  discriminant. 
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